Ignore:
Timestamp:
Dec 16, 2024, 5:30:48 PM (6 days ago)
Author:
jbclement
Message:

PEM:
Addition of the shifting of the soil temperature profile to follow the surface evolution due to ice condensation/sublimation + small cleanings/improvements.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3543 r3553  
    516516== 10/12/2024 == JBC
    517517Using the correct subroutine to output 'ice_porefilling'.
     518
     519== 16/12/2024 == JBC
     520Addition 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  
    216216END SUBROUTINE ini_tsoil_pem
    217217
     218!=======================================================================
     219
     220SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil)
     221
     222use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM
     223use math_mod,      only: solve_steady_heat
     224
     225implicit 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! -------
     233integer,                       intent(in) :: ngrid       ! number of (horizontal) grid-points
     234integer,                       intent(in) :: nsoil       ! number of soil layers
     235integer,                       intent(in) :: nslope      ! number of sub-slopes
     236real, dimension(ngrid,nslope), intent(in) :: zshift_surf ! elevation shift for the surface [m]
     237real, dimension(ngrid,nslope), intent(in) :: zlag        ! newly built lag thickness [m]
     238real, dimension(ngrid,nslope), intent(in) :: tsurf       ! surface temperature [K]
     239! Outputs:
     240! --------
     241real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K]
     242! Local:
     243! ------
     244integer                             :: ig, isoil, islope, ishift, ilag, iz
     245real                                :: a, z, zshift_surfloc
     246real, dimension(ngrid,nsoil,nslope) :: tsoil_old
     247
     248write(*,*)"Shifting soil temperature to surface"
     249tsoil_old = tsoil
     250
     251do 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
     330enddo
     331
     332END SUBROUTINE shift_tsoil2surf
     333
    218334END MODULE compute_soiltemp_mod

  • trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90

    r3498 r3553  
    77!=======================================================================
    88
    9 SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice)
     9SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf)
    1010
    1111use time_evol_mod, only: dt
     12use glaciers_mod,  only: rho_co2ice
    1213
    1314implicit none
    … …  
    1819!
    1920!=======================================================================
    20 !   arguments:
    21 !   ----------
    22 !   INPUT
     21! arguments:
     22! ----------
     23! INPUT
    2324integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes
    2425
    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
     27real, dimension(ngrid,nslope), intent(inout) :: co2_ice     ! Previous and actual density of CO2 ice
     28real, dimension(ngrid,nslope), intent(inout) :: d_co2ice    ! Evolution of perennial ice over one year
     29real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
    2830
    29 !   local:
    30 !   ------
     31! local:
     32! ------
    3133real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice
    3234!=======================================================================
    … …  
    4143end where
    4244
     45zshift_surf = d_co2ice*dt/rho_co2ice
     46
    4347END SUBROUTINE evol_co2_ice
    4448
    4549!=======================================================================
    4650
    47 SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)
     51SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM)
    4852
    49 use time_evol_mod,          only: dt
    50 use comslope_mod,           only: subslope_dist, def_slope_mean
     53use time_evol_mod, only: dt
     54use comslope_mod,  only: subslope_dist, def_slope_mean
     55use glaciers_mod,  only: rho_h2oice
    5156
    5257#ifndef CPP_STD
    … …  
    6368!
    6469!=======================================================================
    65 !   arguments:
    66 !   ----------
    67 !   INPUT
     70! arguments:
     71! ----------
     72! INPUT
    6873integer,                intent(in) :: ngrid, nslope            ! # of grid points, # of subslopes
    6974real, dimension(ngrid), intent(in) :: cell_area                ! Area of each mesh grid (m^2)
    … …  
    7176real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2)
    7277
    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
     79real, dimension(ngrid,nslope), intent(inout) :: h2o_ice     ! Previous and actual density of h2o ice (kg/m^2)
     80real, dimension(ngrid,nslope), intent(inout) :: d_h2oice    ! Evolution of perennial ice over one year (kg/m^2/year)
     81integer,                       intent(inout) :: stopPEM     ! Stopping criterion code
     82real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m]
    7783
    78 !   local:
    79 !   ------
     84! local:
     85! ------
    8086integer                       :: i, islope                                           ! Loop variables
    8187real                          :: 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 done
     88real, dimension(ngrid,nslope) :: new_tend                                            ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done
    8389!=======================================================================
    8490if (ngrid /= 1) then ! Not in 1D
    … …  
    164170endif
    165171
     172zshift_surf = d_h2oice*dt/rho_h2oice
     173
    166174END SUBROUTINE evol_h2o_ice
    167175

  • trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90

    r3532 r3553  
    1010
    1111! 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]
     12real :: inf_h2oice_threshold   ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2]
     13real :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2]
     14real :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2]
     15
     16real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]
     17real, parameter :: rho_h2oice = 920.  ! Density of H2O ice [kg.m-3]
    1518
    1619!=======================================================================

  • trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90

    r3498 r3553  
    88!=======================================================================
    99
     10use glaciers_mod, only: rho_co2ice, rho_h2oice
    1011
    1112implicit none
    … …  
    2324real, parameter :: co2ice_porosity = 0.3        ! Porosity of CO2 ice
    2425real, 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]
    2726real, parameter :: rho_dust = 2500.             ! Density of dust [kg.m-3]
    2827real, parameter :: rho_rock = 3200.             ! Density of rock [kg.m-3]
    … …  
    452451!=======================================================================
    453452! 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
     453SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,zshift_surf,new_lag1,new_lag2,zlag,stopPEM,current1,current2)
     454
     455implicit none
     456
     457!---- Arguments
     458! Inputs
     459real, intent(in) :: d_co2ice, d_h2oice, d_dust
     460
     461! Outputs
    459462type(layering),         intent(inout) :: this
    460463type(stratum), pointer, intent(inout) :: current1, current2
    461464logical,                intent(inout) :: new_str, new_lag1, new_lag2
    462465integer,                intent(inout) :: stopPEM
    463 real,     intent(in) :: d_co2ice, d_h2oice, d_dust
     466real, intent(out) :: zshift_surf, zlag
    464467
    465468!---- Local variables
    … …  
    474477dh_h2oice = d_h2oice/rho_h2oice
    475478dh_dust = d_dust/rho_dust
     479zshift_surf = 0.
     480zlag = 0.
    476481
    477482if (dh_dust < 0.) then ! Dust lifting only
    … …  
    504509        stopPEM = 8
    505510    endif
     511    zshift_surf = dh_dust + h2lift_tot
    506512
    507513else
    … …  
    521527             endif
    522528        endif
     529        zshift_surf = thickness
    523530
    524531!------ Condensation of CO2 ice + H2O ice
    … …  
    536543             endif
    537544        endif
     545        zshift_surf = thickness
    538546
    539547!------ Sublimation of CO2 ice + Condensation of H2O ice
    … …  
    563571                h2subl = h2subl_tot
    564572                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)
    566574            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
    567575                ! Is the considered stratum a dust lag?
    … …  
    575583                h2subl = h_co2ice_old
    576584                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)
    578586                current1 => current1%down
    579587                new_lag1 = .true.
    … …  
    595603             endif
    596604        endif
     605        zshift_surf = zshift_surf + thickness
    597606
    598607!------ Sublimation of CO2 ice + H2O ice
    … …  
    622631                h2subl = h2subl_tot
    623632                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)
    625634            else if (h_co2ice_old < eps) then ! There is nothing to sublimate
    626635                ! Is the considered stratum a dust lag?
    … …  
    634643                h2subl = h_co2ice_old
    635644                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)
    637646                current1 => current1%down
    638647                new_lag1 = .true.
    … …  
    666675                h2subl = h2subl_tot
    667676                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)
    669678            else if (h_h2oice_old < eps) then ! There is nothing to sublimate
    670679                ! Is the considered stratum a dust lag?
    … …  
    678687                h2subl = h_h2oice_old
    679688                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)
    681690                current2 => current2%down
    682691                new_lag2 = .true.
    … …  
    698707!=======================================================================
    699708! 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)
     709SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)
    701710
    702711implicit none
    … …  
    706715type(stratum), pointer, intent(inout) :: str
    707716logical,                intent(inout) :: new_lag
     717real,                   intent(inout) :: zshift_surf, zlag
    708718real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
    709719
    … …  
    719729
    720730str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust
     731zshift_surf = zshift_surf - h2subl/(1. - co2ice_porosity) - hsubl_dust
    721732if (str%thickness < tol) then
    722733    ! Remove the sublimating stratum if there is no ice anymore
    … …  
    748759    endif
    749760endif
     761zlag = zlag + h_lag
    750762
    751763END SUBROUTINE subl_co2ice_layering
    … …  
    753765!=======================================================================
    754766! 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)
     767SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag)
    756768
    757769implicit none
    … …  
    761773type(stratum), pointer, intent(inout) :: str
    762774logical,                intent(inout) :: new_lag
     775real,                   intent(inout) :: zshift_surf, zlag
    763776real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old
    764777
    … …  
    774787
    775788str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust
     789zshift_surf = zshift_surf - h2subl/(1. - h2oice_porosity) - hsubl_dust
    776790if (str%thickness < tol) then
    777791    ! Remove the sublimating stratum if there is no ice anymore
    … …  
    803817    endif
    804818endif
     819zlag = zlag + h_lag
    805820
    806821END SUBROUTINE subl_h2oice_layering

  • trunk/LMDZ.COMMON/libf/evolution/math_mod.F90

    r3532 r3553  
    2828implicit none
    2929
    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!-------
     32integer,             intent(in) :: nz       ! number of layer
     33real, dimension(nz), intent(in) :: z        ! depth layer
     34real, dimension(nz), intent(in) :: y        ! function which needs to be derived
     35real,                intent(in) :: y0, ybot ! boundary conditions
     36! Outputs
     37!--------
     38real, dimension(nz), intent(out) :: dzY ! derivative of y w.r.t depth
     39! Local variables
     40!----------------
    3641integer :: j
    3742real    :: hm, hp, c1, c2, c3
    … …  
    6974implicit none
    7075
    71 integer, intent(in) :: nz
    72 real,    intent(in) :: z(nz), y(nz), y0, yNp1
    73 real, intent(out) :: yp2(nz)
     76! Inputs
     77!-------
     78integer,             intent(in) :: nz
     79real, dimension(nz), intent(in) :: z, y
     80real,                intent(in) :: y0, yNp1
     81! Outputs
     82!--------
     83real, dimension(nz), intent(out) :: yp2
     84! Local variables
     85!----------------
    7486integer :: j
    7587real    :: hm, hp, c1, c2, c3
    … …  
    102114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    103115
    104   implicit none
    105 
    106 integer, intent(in) :: nz, j
    107 real,    intent(in) :: z(nz), y(nz)
     116implicit none
     117
     118! Inputs
     119!-------
     120integer,             intent(in) :: nz, j
     121real, dimension(nz), intent(in) :: z, y
     122! Outputs
     123!--------
    108124real, intent(out) :: dy_zj
     125! Local viariables
     126!-----------------
    109127real :: h1, h2, c1, c2, c3
    110128
    … …  
    135153implicit none
    136154
    137 integer, intent(in) :: nz, i1, i2
    138 real,    intent(in) :: y(nz), z(nz)
     155! Inputs
     156!-------
     157integer,             intent(in) :: nz, i1, i2
     158real, dimension(nz), intent(in) :: y, z
     159! Outputs
     160!--------
    139161real, intent(out) :: integral
    140 integer i
    141 real dz(nz)
     162! Local viariables
     163!-----------------
     164integer             :: i
     165real, dimension(nz) :: dz
    142166
    143167dz(1) = (z(2) - 0.)/2
    … …  
    163187implicit none
    164188
     189! Inputs
     190!-------
    165191real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
    166192real, intent(in) :: z1, z2 ! depth [m]
    167 real, intent(out) :: zr    ! depth at which we have zero
     193! Outputs
     194!--------
     195real, intent(out) :: zr ! depth at which we have zero
    168196
    169197zr = (y1*z2 - y2*z1)/(y1 - y2)
    … …  
    171199END SUBROUTINE findroot
    172200
     201!=======================================================================
     202
     203SUBROUTINE 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
     217implicit none
     218
     219! Inputs
     220!-------
     221integer,                intent(in) :: n
     222real, dimension(n),     intent(in) :: b, d
     223real, dimension(n - 1), intent(in) :: a, c
     224! Outputs
     225!--------
     226real, dimension(n), intent(out) :: x
     227integer,            intent(out) :: error
     228! Local viariables
     229!-----------------
     230integer            :: i
     231real               :: m
     232real, dimension(n) :: cp, dp
     233
     234! Check stability: diagonally dominant condition
     235error = 0
     236if (abs(b(1)) < abs(c(1))) then
     237    error = 1
     238    return
     239endif
     240do i = 2,n - 1
     241    if (abs(b(i)) < abs(a(i - 1)) + abs(c(i))) then
     242        error = 1
     243        return
     244    endif
     245enddo
     246if (abs(b(n)) < abs(a(n - 1))) then
     247    error = 1
     248    return
     249endif
     250
     251! Initialization
     252cp(1) = c(1)/b(1)
     253dp(1) = d(1)/b(1)
     254
     255! Forward sweep
     256do 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
     260enddo
     261m = b(n) - a(n - 1)*cp(n - 1)
     262dp(n) = (d(n) - a(n - 1)*dp(n - 1))/m
     263
     264! Backward substitution
     265x(n) = dp(n)
     266do i = n - 1,1,-1
     267    x(i) = dp(i) - cp(i)*x(i + 1)
     268enddo
     269
     270END SUBROUTINE solve_tridiag
     271
     272!=======================================================================
     273
     274SUBROUTINE 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
     290use abort_pem_mod, only: abort_pem
     291
     292implicit none
     293
     294! Inputs
     295!-------
     296integer,                intent(in) :: n
     297real, dimension(n),     intent(in) :: z, kappa
     298real, dimension(n - 1), intent(in) :: mz, mkappa
     299real,                   intent(in) :: T_left, q_right
     300! Outputs
     301!--------
     302real, dimension(n), intent(out) :: T
     303! Local viariables
     304!-----------------
     305integer                :: i, error
     306real, dimension(n)     :: b, d
     307real, dimension(n - 1) :: a, c
     308
     309! Initialization
     310a = 0.; b = 0.; c = 0.; d = 0.
     311
     312! Left boundary condition (Dirichlet: prescribed temperature)
     313b(1) = 1.
     314d(1) = T_left
     315
     316! Internal points
     317do 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))
     321enddo
     322
     323! Right boundary condition (Neumann: prescribed temperature)
     324a(n - 1) = kappa(n - 1)/(z(n) - z(n - 1))
     325b(n) = -kappa(n)/(z(n) - z(n - 1))
     326d(n) = q_right
     327
     328! Solve the tridiagonal linear system with the Thomas' algorithm
     329call solve_tridiag(a,b,c,d,n,T,error)
     330if (error /= 0) call abort_pem("solve_steady_heat","Unstable solving!",1)
     331
     332END SUBROUTINE solve_steady_heat
     333
    173334end module math_mod

  • trunk/LMDZ.COMMON/libf/evolution/pem.F90

    r3543 r3553  
    4444use constants_marspem_mod,      only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2
    4545use 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, &
     46use comsoil_h_PEM,              only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, mlayer_PEM, &
    4747                                      TI_PEM,               & ! soil thermal inertia
    4848                                      tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
    … …  
    5050use adsorption_mod,             only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
    5151                                      ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays
    52                                       co2_adsorbded_phys, h2o_adsorbded_phys        ! Mass of co2 and h2O adsorbded
     52                                      co2_adsorbed_phys, h2o_adsorbed_phys        ! Mass of co2 and h2O adsorbed
    5353use time_evol_mod,              only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini
    5454use orbit_param_criterion_mod,  only: orbit_param_criterion
    … …  
    6565use pemetat0_mod,               only: pemetat0
    6666use read_data_PCM_mod,          only: read_data_PCM
    67 use recomp_tend_co2_slope_mod,  only: recomp_tend_co2_slope
    68 use compute_soiltemp_mod,       only: compute_tsoil_pem
     67use recomp_tend_co2_slope_mod,  only: recomp_tend_co2
     68use compute_soiltemp_mod,       only: compute_tsoil_pem, shift_tsoil2surf
    6969use writediagpem_mod,           only: writediagpem, writediagsoilpem
    7070use co2condens_mod,             only: CO2cond_ps
    … …  
    217217real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
    218218real, dimension(:,:),     allocatable :: Tsurfavg_before_saved            ! Surface temperature saved from previous time step [K]
    219 real, dimension(:),       allocatable :: delta_co2_adsorbded              ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
    220 real, dimension(:),       allocatable :: delta_h2o_adsorbded              ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
    221 real                                  :: totmassco2_adsorbded             ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
    222 real                                  :: totmassh2o_adsorbded             ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
     219real, dimension(:),       allocatable :: delta_co2_adsorbed              ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
     220real, dimension(:),       allocatable :: delta_h2o_adsorbed              ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
     221real                                  :: totmassco2_adsorbed             ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
     222real                                  :: totmassh2o_adsorbed             ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
    223223logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
    224224logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
    … …  
    228228real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
    229229real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
     230real, dimension(:,:),     allocatable :: zshift_surf                      ! Elevation shift for the surface [m]
     231real, dimension(:,:),     allocatable :: zlag                             ! Newly built lag thickness [m]
    230232
    231233! Some variables for the PEM run
    … …  
    541543allocate(vmr_co2_PCM(ngrid,timelen))
    542544allocate(ps_timeseries(ngrid,timelen))
    543 allocate(min_co2_ice(ngrid,nslope,2))
    544 allocate(min_h2o_ice(ngrid,nslope,2))
     545allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2))
    545546allocate(tsurf_avg_yr1(ngrid,nslope))
    546547allocate(tsurf_avg(ngrid,nslope))
    547548allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen))
    548549allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen))
    549 allocate(q_co2_PEM_phys(ngrid,timelen))
    550 allocate(q_h2o_PEM_phys(ngrid,timelen))
     550allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope))
     551allocate(q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen))
    551552allocate(watersurf_density_avg(ngrid,nslope))
    552553allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen))
    … …  
    554555allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    555556allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen))
    556 allocate(delta_co2_adsorbded(ngrid))
     557allocate(delta_co2_adsorbed(ngrid))
    557558allocate(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))
     559allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope))
     560allocate(delta_h2o_icetablesublim(ngrid),delta_h2o_adsorbed(ngrid))
    562561allocate(vmr_co2_PEM_phys(ngrid,timelen))
    563562
    … …  
    649648              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
    650649              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_adsorbded_phys,         &
    652               delta_co2_adsorbded,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)
    653652
    654653! We save the places where h2o ice is sublimating
    … …  
    680679
    681680if (adsorption_pem) then
    682     totmassco2_adsorbded = 0.
    683     totmassh2o_adsorbded = 0.
     681    totmassco2_adsorbed = 0.
     682    totmassh2o_adsorbed = 0.
    684683    do ig = 1,ngrid
    685684        do islope = 1,nslope
    686685            do l = 1,nsoilmx_PEM - 1
    687686                if (l==1) then
    688                    totmassco2_adsorbded = 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))* &
    689688                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    690                    totmassh2o_adsorbded = 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))* &
    691690                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    692691                else
    693                    totmassco2_adsorbded = 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))* &
    694693                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    695                    totmassh2o_adsorbded = 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))* &
    696695                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    697696                endif
    … …  
    699698        enddo
    700699    enddo
    701     write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbded
    702     write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbded
     700    write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbed
     701    write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed
    703702endif ! adsorption
    704703deallocate(tsurf_avg_yr1)
    … …  
    750749    if (adsorption_pem) then
    751750        do i = 1,ngrid
    752             global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
     751            global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbed(i)/Total_surface
    753752        enddo
    754753    endif
    … …  
    819818!    II_b Evolution of ice
    820819!------------------------
    821     call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,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)
    823822    if (layering_algo) then
    824823        do islope = 1,nslope
    825824            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)
    827826            enddo
    828827        enddo
     828    else
     829        zlag = 0.
    829830    endif
    830831
    … …  
    887888                tsurf_avg(ig,islope) = ave/timelen
    888889            endif
     890            !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBLIMATION??? !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    889891        enddo
    890892    enddo
    … …  
    901903
    902904    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
    905909        write(*,*)"Updating soil temperature"
    906910        allocate(Tsoil_locslope(ngrid,nsoilmx_PEM))
    … …  
    926930        deallocate(Tsoil_locslope)
    927931
    928 ! II_d.3 Update the ice table
     932! II_d.4 Update the ice table
    929933        if (icetable_equilibrium) then
    930934            write(*,*) "Compute ice table (equilibrium method)"
    … …  
    947951        endif
    948952
    949 ! II_d.4 Update the soil thermal properties
     953! II_d.5 Update the soil thermal properties
    950954        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)
    951955
    952 ! II_d.5 Update the mass of the regolith adsorbed
     956! II_d.6 Update the mass of the regolith adsorbed
    953957        if (adsorption_pem) then
    954958            call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice,                           &
    955959                                     h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, &
    956                                      h2o_adsorbded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)
    957 
    958             totmassco2_adsorbded = 0.
    959             totmassh2o_adsorbded = 0.
     960                                     h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed)
     961
     962            totmassco2_adsorbed = 0.
     963            totmassh2o_adsorbed = 0.
    960964            do ig = 1,ngrid
    961965                do islope = 1,nslope
    962966                    do l = 1,nsoilmx_PEM
    963967                        if (l == 1) then
    964                             totmassco2_adsorbded = 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))* &
    965969                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    966                             totmassh2o_adsorbded = 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))* &
    967971                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    968972                        else
    969                             totmassco2_adsorbded = 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))* &
    970974                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    971                             totmassh2o_adsorbded = 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))* &
    972976                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
    973977                        endif
    … …  
    975979                enddo
    976980            enddo
    977             write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbded
    978             write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbded
     981            write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbed
     982            write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbed
    979983        endif
    980984    endif !soil_pem
    … …  
    10051009            if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope))
    10061010            if (adsorption_pem) then
    1007                 call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbded_phys(:,:,islope))
    1008                 call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbded_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))
    10091013            endif
    10101014        endif
    … …  
    10151019!    II_f Update the tendencies
    10161020!------------------------
    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)
    10181022
    10191023!------------------------
    … …  
    10211025!    II_g Checking the stopping criterion
    10221026!------------------------
    1023 
    10241027    write(*,*) "Checking the stopping criteria..."
    10251028    call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid)
    … …  
    12101213call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
    12111214call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
    1212              co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
     1215             co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif)
    12131216write(*,*) "restartpem.nc has been written"
    12141217
    … …  
    12311234deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
    12321235deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
    1233 deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
    1234 deallocate(icetable_thickness_old,ice_porefilling_old)
     1236deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,vmr_co2_PEM_phys,delta_h2o_icetablesublim)
     1237deallocate(icetable_thickness_old,ice_porefilling_old,zshift_surf,zlag)
    12351238deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
    12361239!----------------------------- END OUTPUT ------------------------------

  • trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90

    r3498 r3553  
    6565write(*,*) 'recomp_orb_param, Old ecc. =', e_elips
    6666write(*,*) 'recomp_orb_param, Old Lsp  =', Lsp
     67write(*,*) 'recomp_orb_param, New year =', Year
    6768
    6869found_year = .false.
    … …  
    116117#endif
    117118
    118 write(*,*) 'recomp_orb_param, New year =', Year
    119119write(*,*) 'recomp_orb_param, New obl. =', obliquit
    120120write(*,*) 'recomp_orb_param, New ecc. =', e_elips

  • trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90

    r3532 r3553  
    77!=======================================================================
    88
    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)
     9SUBROUTINE 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)
    1111
    1212use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol
    … …  
    2424!   INPUT
    2525integer,                        intent(in) :: timelen, ngrid, nslope
    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_phys_ini ! physical point field: Evolution of perennial ice over one year
    32 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)
     26real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM          ! physical point field: Volume mixing ratio of co2 in the first layer
     27real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM          ! physical point field: Volume mixing ratio of co2 in the first layer
     28real, dimension(ngrid,timelen), intent(in) :: ps_PCM_2             ! physical point field: Surface pressure in the PCM
     29real,                           intent(in) :: global_avg_press_PCM ! global averaged pressure at previous timestep
     30real,                           intent(in) :: global_avg_press_new ! global averaged pressure at current timestep
     31real, dimension(ngrid,nslope),  intent(in) :: d_co2ice_ini        ! physical point field: Evolution of perennial ice over one year
     32real, dimension(ngrid,nslope),  intent(in) :: co2ice               ! CO2 ice per mesh and sub-grid slope (kg/m^2)
     33real, dimension(ngrid,nslope),  intent(in) :: emissivity           ! Emissivity per mesh and sub-grid slope(1)
    3434!   OUTPUT
    35 real, dimension(ngrid,nslope), intent(inout) ::  d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
     35real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year
    3636
    3737!   local:
    … …  
    4545do i = 1,ngrid
    4646    do islope = 1,nslope
    47         coef = sols_per_my*sec_per_sol*emissivity_slope(i,islope)*sigmaB/Lco2
     47        coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2
    4848        ave = 0.
    49         if (co2ice_slope(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
    50             do t=1,timelen
     49        if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then
     50            do t = 1,timelen
    5151                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.)))**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.)))**4
    5353            enddo
    54             if (ave < 1e-4) ave = 0.
    55             d_co2ice_phys(i,islope) = d_co2ice_phys_ini(i,islope) - coef*ave/timelen
     54            if (ave < 1.e-4) ave = 0.
     55            d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*ave/timelen
    5656        endif
    5757    enddo
    5858enddo
    5959
    60 END SUBROUTINE recomp_tend_co2_slope
     60END SUBROUTINE recomp_tend_co2
     61!=======================================================================
     62
     63SUBROUTINE recomp_tend_h2o(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp)
     64
     65implicit none
     66
     67!=======================================================================
     68!
     69!  Routine that compute the evolution of the tendencie for h2o ice
     70!
     71!=======================================================================
     72
     73!   arguments:
     74!   ----------
     75!   INPUT
     76integer,            intent(in) :: timelen, ngrid, nslope
     77real, dimension(:), intent(in) :: PCM_temp, PEM_temp
     78
     79!   OUTPUT
     80real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year
     81
     82!   local:
     83!   ------
     84real :: coef, ave, Rz_old, Rz_new, R_dec, soil_psv_old, soil_psv_new, hum_dec, h2oice_depth_new, h2oice_depth_old
     85
     86write(*,*) "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
     97END SUBROUTINE recomp_tend_h2o
    6198
    6299END MODULE recomp_tend_co2_slope_mod
Note: See TracChangeset for help on using the changeset viewer.