Changeset 2961


Ignore:
Timestamp:
May 9, 2023, 12:18:41 PM (20 months ago)
Author:
llange
Message:

PEM

  • Move math functions (findroot, deriv, etc) in a specific module
  • Ice table is no longer infinite, the bottom of the ice table (zend) is set such at d(rho_vap) dz (zend) = 0. Future commit will transfer the loss of ice table to perenial glaciers

LL

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
1 added
1 deleted
8 edited

Legend:

Unmodified
Added
Removed

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

    r2918 r2961  
    1919  USE adsorption_mod,only: adsorption_pem
    2020  USE co2glaciers_mod, only: co2glaciersflow
     21  use ice_table_mod, only: icetable_equilibrium, icetable_dynamic
    2122
    2223  CHARACTER(len=20),parameter :: modname ='conf_pem'
    … …  
    8485  print*, 'Thermal properties of the regolith vary with pressure ?', reg_thprop_dependp
    8586
     87   icetable_equilibrium = .true.
     88   CALL getin('icetable_equilibrium',icetable_equilibrium)
     89   print*, 'Do we compute the ice table at equilibrium?', icetable_equilibrium
     90
     91   icetable_dynamic = .false.
     92   CALL getin('icetable_dynamic',icetable_dynamic)
     93   print*, 'Do we compute the ice table with the dynamic method?', icetable_dynamic
     94
     95  if ((not(soil_pem)).and.((icetable_equilibrium).or.(icetable_dynamic))) then
     96       print*,'Ice table  must be used when soil_pem = T'
     97       call abort_physic(modname,"Ice table  must be used when soil_pem = T",1)
     98  endif
    8699
    87100  if ((not(soil_pem)).and.adsorption_pem) then
    … …  
    90103  endif
    91104 
    92   if ((not(soil_pem)).and.fluxgeo.gt.0.) then
     105  if ((not(soil_pem)).and.(fluxgeo.gt.0.)) then
    93106       print*,'Soil is not activated but Flux Geo > 0.'
    94107       call abort_physic(modname,"Soil is not activated but Flux Geo > 0.",1)

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

    r2945 r2961  
    2828      REAL,PARAMETER :: TI_bedrock = 2300.            ! Thermal inertia of Bedrock following Wood 2009 [SI]
    2929
     30!     Porosity of the soil
     31      REAL,PARAMETER :: porosity = 0.4                ! porosity of the martian soil, correspond to the value for a random loose packing of monodiperse sphere (Scott, 1960)
     32
    3033!     Stefan Boltzmann constant
    3134      REAL,PARAMETER :: sigmaB=5.678E-8
    … …  
    3538
    3639
     40     
    3741    END MODULE constants_marspem_mod
    3842

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

    r2902 r2961  
    33  implicit none
    44
     5!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     6!!!
     7!!! Purpose: Ice table (pore-filling) variables and methods to compute it (dynamic and static)
     8!!! Author:  LL, 02/2023
     9!!!
     10!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     11
     12  LOGICAL,SAVE :: icetable_equilibrium                    ! Boolean to say if the PEM needs to recompute the icetable depth when at  equilibrium
     13  LOGICAL,SAVE :: icetable_dynamic                        ! Boolean to say if the PEM needs to recompute the icetable depth (dynamic method)
     14  real,save,allocatable  :: porefillingice_depth(:,:)    ! ngrid x nslope: Depth of the ice table [m]
     15  real,save,allocatable  :: porefillingice_thickness(:,:)    ! ngrid x nslope: Thickness of the ice table [m]
     16
    517  contains
    618
    719
    8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    9 !!!
    10 !!! Purpose: Compute the ice table in two ways: dynamic and at equilibrium
    11 !!! Author:  LL, 02/2023
    12 !!!
    13 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    14 
    15 
    16 
    17    SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,ice_table)
     20  subroutine ini_ice_table_porefilling(ngrid,nslope)
     21 
     22    implicit none
     23    integer,intent(in) :: ngrid ! number of atmospheric columns
     24    integer,intent(in) :: nslope ! number of slope within a mesh
     25
     26    allocate(porefillingice_depth(ngrid,nslope))
     27    allocate(porefillingice_thickness(ngrid,nslope))
     28
     29  end subroutine ini_ice_table_porefilling
     30
     31  subroutine end_ice_table_porefilling
     32
     33    implicit none
     34    if (allocated(porefillingice_depth)) deallocate(porefillingice_depth)
     35    if (allocated(porefillingice_thickness)) deallocate(porefillingice_thickness)
     36
     37  end subroutine end_ice_table_porefilling
     38
     39
     40
     41
     42   SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
    1843!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1944!!!
    20 !!! Purpose: Compute the ice table depth knowing the yearly average water
     45!!! Purpose: Compute the ice table depth (pore-filling) knowing the yearly average water
    2146!!! density at the surface and at depth.
    2247!!! Computations are made following the methods in Schorgofer et al., 2005
    … …  
    2651!!!
    2752!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    28 
     53    use math_mod,only: findroot
    2954#ifndef CPP_STD
    30     USE comsoil_h_PEM, only: mlayer_PEM                             ! Depth of the vertical grid
     55    USE comsoil_h_PEM, only: mlayer_PEM,layer_PEM                  ! Depth of the vertical grid
     56    USE soil_thermalproperties_mod, only: ice_thermal_properties
    3157    implicit none
    32 
     58!   inputs
     59! -----------
    3360    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM
    3461    logical,intent(in) :: watercaptag(ngrid)                       ! Boolean to check the presence of a perennial glacier
    3562    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
    3663    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
    37     real,intent(inout) :: ice_table(ngrid,nslope)                  ! ice table depth [m]
     64    real,intent(in) :: regolith_inertia(ngrid,nslope)              ! Thermal inertia of the regolith layer [SI]
     65!   Ouputs
     66! -----------
     67    real,intent(out) :: ice_table_beg(ngrid,nslope)               ! ice table depth [m]
     68    real,intent(out) :: ice_table_thickness(ngrid,nslope)         ! ice table thickness [m]
     69!   Local
     70! -----------
    3871    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root
    39     integer ig, islope,isoil                                       ! loop variables
     72    integer ig, islope,isoil,isoilend                              ! loop variables
    4073    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]
    41 
    42 
     74    real  :: ice_table_end                                         ! depth of the end of the ice table  [m]
     75    real :: previous_icetable_depth(ngrid,nslope)                  ! Ice table computed at previous ice depth [m]
     76    real :: stretch                                                ! stretch factor to improve the convergence of the ice table
     77    real :: wice_inertia                                           ! Water Ice thermal Inertia [USI]   
     78!   Code
     79! -----------
     80
     81     previous_icetable_depth(:,:) = ice_table_beg(:,:)
    4382     do ig = 1,ngrid
    4483      if(watercaptag(ig)) then
    45          ice_table(ig,:) = 0.
     84         ice_table_beg(ig,:) = 0.
     85         ice_table_thickness(ig,:) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
    4686      else
    4787        do islope = 1,nslope
    48            ice_table(ig,islope) = -1.
    49 
     88           ice_table_beg(ig,islope) = -1.
     89           ice_table_thickness(ig,islope) = 0.
    5090         do isoil = 1,nsoil_PEM
    5191           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
    52 
    5392         enddo
    5493         if(diff_rho(1) > 0) then                                   ! ice is at the surface
    55            ice_table(ig,islope) = 0.
     94           ice_table_beg(ig,islope) = 0.
     95               do isoilend = 2,nsoil_PEM-1
     96                 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
     97                  call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
     98                  ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
     99                  exit
     100                 endif
     101               enddo
    56102         else
    57103           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
    58104             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
    59                call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table(ig,islope))
    60                exit
     105               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table_beg(ig,islope))
     106               ! Now let's find the end of the ice table
     107               ice_table_thickness(ig,islope) = layer_PEM(nsoil_PEM) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     108               do isoilend = isoil+1,nsoil_PEM-1
     109                 if((diff_rho(isoilend).gt.0).and.(diff_rho(isoilend+1).lt.0.)) then
     110                  call findroot(diff_rho(isoilend),diff_rho(isoilend+1),mlayer_PEM(isoilend),mlayer_PEM(isoilend+1),ice_table_end)
     111                  ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
     112                  exit
     113                 endif
     114               enddo
    61115             endif !diff_rho(z) <0 & diff_rho(z+1) > 0
    62116           enddo
    … …  
    65119       endif ! watercaptag
    66120      enddo
     121
     122! Small trick to speed up the convergence, Oded's idea.
     123      do islope = 1,nslope
     124        do ig = 1,ngrid
     125         if((ice_table_beg(ig,islope).gt.previous_icetable_depth(ig,islope)).and.(previous_icetable_depth(ig,islope).ge.0)) then
     126            call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
     127            stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
     128
     129            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
     130            previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
     131            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope)+(ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
     132         endif
     133        enddo
     134      enddo
     135
    67136!=======================================================================
    68137      RETURN
    … …  
    70139      END
    71140
     141
    72142   SUBROUTINE find_layering_icetable(porefill,psat_soil,psat_surf,pwat_surf,psat_bottom,B,index_IS,depth_filling,index_filling,index_geothermal,depth_geothermal,dz_etadz_rho)
    73143!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    79149!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    80150use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM
     151use math_mod, only: deriv1,deriv1_onesided,colint,findroot,deriv2_simple
    81152implicit none
    82153! inputs
    … …  
    232303
    233304
    234 
    235 SUBROUTINE findroot(y1,y2,z1,z2,zr)
    236         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    237         !!!
    238         !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
    239         !!!
    240         !!! Author: LL
    241         !!!
    242         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    243 
    244         implicit none
    245         real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
    246         real,intent(in) :: z1,z2 ! depth [m]
    247         real,intent(out) :: zr ! depth at which we have zero
    248         zr = (y1*z2 - y2*z1)/(y1-y2)
    249         RETURN
    250         end
    251 
    252305SUBROUTINE constriction(porefill,eta)
    253306
    … …  
    272325        if (porefill.le.1.) eta = 0.
    273326        return
    274         end
    275 
    276 
    277 
    278 subroutine deriv1(z,nz,y,y0,ybot,dzY)
    279 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    280 !!!
    281 !!! Purpose: Compute the first derivative of a function y(z) on an irregular grid
    282 !!!           
    283 !!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
    284 !!!
    285 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    286         implicit none
    287 
    288 ! first derivative of a function y(z) on irregular grid
    289 ! upper boundary conditions: y(0)=y0
    290 ! lower boundary condition.: yp = ybottom
    291   integer, intent(IN) :: nz ! number of layer
    292   real, intent(IN) :: z(nz) ! depth layer
    293   real, intent(IN) :: y(nz) ! function which needs to be derived
    294   real, intent(IN) :: y0,ybot ! boundary conditions
    295   real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
    296 ! local
    297   integer :: j
    298   real :: hm,hp,c1,c2,c3
    299 
    300   hp = z(2)-z(1)
    301   c1 = z(1)/(hp*z(2))
    302   c2 = 1/z(1) - 1/(z(2)-z(1))
    303   c3 = -hp/(z(1)*z(2))
    304   dzY(1) = c1*y(2) + c2*y(1) + c3*y0
    305   do j=2,nz-1
    306      hp = z(j+1)-z(j)
    307      hm = z(j)-z(j-1)
    308      c1 = +hm/(hp*(z(j+1)-z(j-1)))
    309      c2 = 1/hm - 1/hp
    310      c3 = -hp/(hm*(z(j+1)-z(j-1)))
    311      dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
    312   enddo
    313   dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
    314  return
    315 end subroutine deriv1
    316 
    317 
    318 
    319 subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
    320 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    321 !!!
    322 !!! Purpose: Compute the second derivative of a function y(z) on an irregular grid
    323 !!!           
    324 !!! Author: N.S (raw copy/paste from MSIM)
    325 !!!
    326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    327 
    328   ! second derivative y_zz on irregular grid
    329   ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
    330   implicit none
    331   integer, intent(IN) :: nz
    332   real, intent(IN) :: z(nz),y(nz),y0,yNp1
    333   real, intent(OUT) :: yp2(nz)
    334   integer j
    335   real hm,hp,c1,c2,c3
    336 
    337   c1 = +2./((z(2)-z(1))*z(2))
    338   c2 = -2./((z(2)-z(1))*z(1))
    339   c3 = +2./(z(1)*z(2))
    340   yp2(1) = c1*y(2) + c2*y(1) + c3*y0
    341 
    342   do j=2,nz-1
    343      hp = z(j+1)-z(j)
    344      hm = z(j)-z(j-1)
    345      c1 = +2./(hp*(z(j+1)-z(j-1)))
    346      c2 = -2./(hp*hm)
    347      c3 = +2./(hm*(z(j+1)-z(j-1)))
    348      yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
    349   enddo
    350   yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
    351   return
    352 end subroutine deriv2_simple
    353 
    354 
    355 
    356 subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
    357 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    358 !!!
    359 !!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
    360 !!!           
    361 !!! Author: N.S (raw copy/paste from MSIM)
    362 !!!
    363 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    364 
    365   implicit none
    366   integer, intent(IN) :: nz,j
    367   real, intent(IN) :: z(nz),y(nz)
    368   real, intent(out) :: dy_zj
    369   real h1,h2,c1,c2,c3
    370   if (j<1 .or. j>nz-2) then
    371      dy_zj = -1.
    372   else
    373      h1 = z(j+1)-z(j)
    374      h2 = z(j+2)-z(j+1)
    375      c1 = -(2*h1+h2)/(h1*(h1+h2))
    376      c2 =  (h1+h2)/(h1*h2)
    377      c3 = -h1/(h2*(h1+h2))
    378      dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
    379   endif
    380   return
    381 end subroutine deriv1_onesided
    382 
    383 
    384 subroutine  colint(y,z,nz,i1,i2,integral)
    385 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    386 !!!
    387 !!! Purpose:  Column integrates y on irregular grid
    388 !!!           
    389 !!! Author: N.S (raw copy/paste from MSIM)
    390 !!!
    391 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    392  
    393   implicit none
    394   integer, intent(IN) :: nz, i1, i2
    395   real, intent(IN) :: y(nz), z(nz)
    396   real,intent(out) :: integral
    397   integer i
    398   real dz(nz)
    399  
    400   dz(1) = (z(2)-0.)/2
    401   do i=2,nz-1
    402      dz(i) = (z(i+1)-z(i-1))/2.
    403   enddo
    404   dz(nz) = z(nz)-z(nz-1)
    405   integral = sum(y(i1:i2)*dz(i1:i2))
    406 end subroutine colint
    407 
    408 
    409 
    410 
    411 
    412 
    413         end module
     327        end subroutine
     328
     329
     330end module

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

    r2944 r2961  
    104104      use orbit_param_criterion_mod, only : orbit_param_criterion
    105105      use recomp_orb_param_mod, only: recomp_orb_param
    106       use ice_table_mod, only: computeice_table_equilibrium
     106      use ice_table_mod, only: porefillingice_depth,porefillingice_thickness,&
     107                               end_ice_table_porefilling,ini_ice_table_porefilling, &
     108                               computeice_table_equilibrium
     109      use soil_thermalproperties_mod, only: update_soil_thermalproperties
    107110
    108111
    … …  
    223226     REAL, ALLOCATABLE :: tsurf_ave_yr1(:,:) ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
    224227     REAL, ALLOCATABLE :: inertiesoil(:,:)    !Physic x Depth  Thermal inertia of the mesh for restart [SI]
    225 
    226228     REAL, ALLOCATABLE :: TI_GCM(:,:,:) ! Same but for the start
    227 
    228      REAL,ALLOCATABLE  :: ice_depth(:,:)      ! Physic x SLope: Ice table depth [m]
    229229     REAL,ALLOCATABLE  :: TI_locslope(:,:)    ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
    230230     REAL,ALLOCATABLE  :: Tsoil_locslope(:,:) ! Physic x Soil: intermediate when computing Tsoil [K]
    … …  
    535535!---------------------------- Initialisation of the PEM soil and values ---------------------
    536536
    537       call end_comsoil_h_PEM
    538       call ini_comsoil_h_PEM(ngrid,nslope)
    539       call end_adsorption_h_PEM
    540       call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
    541 
    542       allocate(ice_depth(ngrid,nslope))
    543       ice_depth(:,:) = 0.
     537   call end_comsoil_h_PEM
     538   call ini_comsoil_h_PEM(ngrid,nslope)
     539   call end_adsorption_h_PEM
     540   call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
     541   call end_ice_table_porefilling
     542   call ini_ice_table_porefilling(ngrid,nslope)
    544543
    545544   if(soil_pem) then
    … …  
    673672!---------------------------- Read the PEMstart ---------------------
    674673
    675       call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,ice_depth, &
     674      call pemetat0("startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_PEM_yr1,tsoil_PEM,porefillingice_depth,porefillingice_thickness, &
    676675      tsurf_ave_yr1, tsurf_ave, q_co2_PEM_phys, q_h2o_PEM_phys,ps_timeseries,tsoil_phys_PEM_timeseries,&
    677676      tendencies_h2o_ice,tendencies_co2_ice,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM,&
    678677      watersurf_density_ave,watersoil_density_PEM_ave, &
    679678      co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,water_reservoir)
    680 
    681     do ig = 1,ngrid
    682       do islope = 1,nslope
     679 
     680      do ig = 1,ngrid
     681        do islope = 1,nslope
    683682          qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig)*cos(pi*def_slope_mean(islope)/180.)
     683        enddo
    684684      enddo
    685     enddo
    686685
    687686    if(adsorption_pem) then
    … …  
    987986
    988987! II_d.3 Update the ice table
    989      call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,ice_depth)
    990        
     988     call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
    991989     print *, "Update soil propreties"
    992990
    993991! II_d.4 Update the soil thermal properties
    994       call update_soil(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
    995         ice_depth,TI_PEM)
     992      call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_new, &
     993        porefillingice_depth,porefillingice_thickness,TI_PEM)
    996994
    997995! II_d.5 Update the mass of the regolith adsorbded
    … …  
    13341332
    13351333        call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , &
    1336                       tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
     1334                      tsoil_PEM, TI_PEM, porefillingice_depth,porefillingice_thickness, &
     1335                      co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir)
    13371336
    13381337        call info_run_PEM(year_iter, criterion_stop)

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

    r2945 r2961  
    1    SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
     1   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, ice_table_thickness, &
    22                      tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
    33                      global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
    … …  
    1212   USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, &
    1313                                   TI_breccia,TI_bedrock
     14   use soil_thermalproperties_mod, only: update_soil_thermalproperties
    1415#ifndef CPP_STD   
    1516   use surfdat_h, only: watercaptag
    … …  
    4344  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)             ! soil (mid-layer) temperature [K]
    4445  real,intent(inout) :: ice_table(ngrid,nslope)                       ! Ice table depth [m]
     46  real,intent(inout) :: ice_table_thickness(ngrid,nslope)             ! Ice table thickness [m]
    4547  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen)    ! instantaneous soil (mid-layer) temperature [K]
    4648  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed [kg/m^2]
    … …  
    260262      write(*,*)'PEM settings: failed loading <ice_table>'
    261263      write(*,*)'will reconstruct the values of the ice table given the current state'
    262      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
    263      call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
     264     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, TI_PEM(:,1,:),ice_table,ice_table_thickness)
     265     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    264266     do islope = 1,nslope
    265267     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
    … …  
    453455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    454456!c) Ice table   
    455      call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave, ice_table)
    456      call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,TI_PEM)
     457     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
     458     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
    457459     do islope = 1,nslope
    458460     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))

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

    r2895 r2961  
    7171
    7272subroutine pemdem1(filename,year_iter,nsoil_PEM,ngrid,nslope, &
    73                     tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table, &
     73                    tsoil_slope_PEM,inertiesoil_slope_PEM,ice_table_depth,ice_table_thickness, &
    7474                    m_co2_regolith,m_h2o_regolith,water_reservoir)
    7575
    … …  
    9393  real,intent(IN) :: tsoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    9494  real,intent(IN) :: inertiesoil_slope_PEM(ngrid,nsoil_PEM,nslope) !under mesh bining according to slope
    95   real,intent(IN) :: ice_table(ngrid,nslope) !under mesh bining according to slope
     95  real,intent(IN) :: ice_table_depth(ngrid,nslope) !under mesh bining according to slope
     96  real,intent(IN) :: ice_table_thickness(ngrid,nslope) !under mesh bining according to slope
    9697  integer,intent(IN) :: nslope
    9798  real,intent(in) :: m_co2_regolith(ngrid,nsoil_PEM,nslope)
    … …  
    136137ENDDO
    137138
    138   call put_field("ice_table","Depth of ice table", &
    139                  ice_table,year_tot)
     139  call put_field("ice_table_depth","Depth of ice table", &
     140                 ice_table_depth,year_tot)
     141
     142  call put_field("ice_table_thickness","Depth of ice table", &
     143                 ice_table_thickness,year_tot)
    140144
    141145  call put_field("inertiedat_PEM","Thermal inertie of PEM ", &

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

    r2895 r2961  
    3333
    3434! 0. Initialisations and preprocessing step
    35 
    3635   
    3736! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities

  • trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM.F

    r2895 r2961  
    9292! 3. Index for subsurface layering
    9393! ------------------
    94        write(*,*) 'depth=',depth_breccia,depth_bedrock
    9594      index_breccia= 1
    9695      do iloop = 1,nsoil_PEM-1
Note: See TracChangeset for help on using the changeset viewer.