Index: /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 2901)
+++ /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 2902)
@@ -28,5 +28,5 @@
 
 #ifndef CPP_STD
-    USE comsoil_h_PEM, only: layer_PEM                             ! Depth of the vertical grid 
+    USE comsoil_h_PEM, only: mlayer_PEM                             ! Depth of the vertical grid 
     implicit none 
 
@@ -46,5 +46,5 @@
       else
         do islope = 1,nslope
-           ice_table(ig,islope) = -10.
+           ice_table(ig,islope) = -1.
 
          do isoil = 1,nsoil_PEM
@@ -57,7 +57,5 @@
            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
              if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
-               z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
-               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
-               ice_table(ig,islope) = -z2/z1
+               call findroot(diff_rho(isoil),diff_rho(isoil+1),mlayer_PEM(isoil),mlayer_PEM(isoil+1),ice_table(ig,islope))
                exit
              endif !diff_rho(z) <0 & diff_rho(z+1) > 0
@@ -72,5 +70,344 @@
       END
 
-
-
-end module
+   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)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute layering between dry soil, pore filling ice, and ice sheet based on Schorgofer, Icarus (2010). Adapted from NS MSIM 
+!!! 
+!!! Author: LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+use comsoil_h_PEM,only: nsoilmx_PEM,mlayer_PEM 
+implicit none
+! inputs
+! ------
+
+real,intent(in) :: porefill(nsoilmx_PEM) ! Fraction of pore space filled with ice [Unitless] 0 <= f <= 1 for pore ice
+real,intent(in) :: psat_soil(nsoilmx_PEM) ! Soil water pressure at saturation, yearly averaged [Pa]
+real,intent(in) :: psat_surf ! surface water pressure at saturation, yearly averaged [Pa]
+real,intent(in) :: pwat_surf ! Water vapor pressure  at the surface, not necesseraly at saturation, yearly averaged [Pa]
+real,intent(in) :: psat_bottom ! Boundary conditions for soil vapor pressure [Pa]
+real,intent(in) :: B ! constant (Eq. 8 from  Schorgofer, Icarus (2010).) 
+integer, intent(in) :: index_IS ! index of the soil layer where the ice sheet begins [1]
+real, intent(inout) :: depth_filling ! depth where pore filling begins [m]
+
+! outputs
+! -------
+integer,intent(out) :: index_filling ! index where the pore filling begins [1]
+integer, intent(out) :: index_geothermal ! index where the ice table stops [1]
+real, intent(out) :: depth_geothermal ! depth where the ice table stops [m]
+real, intent(out) :: dz_etadz_rho(nsoilmx_PEM) ! \partial z(eta \partial z rho), eta is the constriction, used later for pore filling increase
+
+! local
+
+real :: eta(nsoilmx_PEM)  ! constriction
+integer :: ilay ! index for loop
+real :: old_depth_filling ! depth_filling saved
+real :: dz_psat(nsoilmx_PEM) ! first derivative of the vapor pressure at saturationn
+integer :: index_tmp ! for loop 
+real :: Jdry ! flux trought the dry layer
+real :: Jsat ! flux trought the ice layer
+real :: Jdry_prevlay,Jsat_prevlay ! same but for the previous ice layer
+integer :: index_firstice ! first index where ice appears (i.e., f > 0)
+real :: dz_eta(nsoilmx_PEM) ! \partial z \eta
+real :: dz_eta_low ! same but evaluated at the interface for ice
+real :: dzz_psat(nsoilmx_PEM) ! \partial \partial psat
+real :: massfillabove,massfillafter ! h2O mass above and after index_geothermal
+
+! constant
+real :: pvap2rho = 18.e-3/8.314
+!
+
+
+
+
+  
+! 0. Compute constriction over the layer
+! Within the ice sheet, constriction is set to 0. Elsewhere, constriction =  (1-porefilling)**2
+if (index_IS.lt.0) then
+   index_tmp = nsoilmx_PEM
+   do ilay = 1,nsoilmx_PEM
+       call constriction(porefill(ilay),eta(ilay))
+   enddo
+else
+   index_tmp = index_IS
+   do ilay = 1,index_IS-1
+       call constriction(porefill(ilay),eta(ilay))
+   enddo
+   do ilay = index_IS,nsoilmx_PEM
+       eta(ilay) = 0.
+   enddo
+endif
+
+! 1. Depth at which pore filling occurs. We solve Eq. 9 from  Schorgofer, Icarus  (2010)    
+
+old_depth_filling = depth_filling
+
+call deriv1(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dz_psat)  
+
+do ilay = 1,index_tmp
+  Jdry = (psat_soil(ilay) - pwat_surf)/mlayer_PEM(ilay) ! left member of Eq. 9 from Schorgofer, Icarus  (2010)
+  Jsat =  eta(ilay)*dz_psat(ilay) !right member of Eq. 9 from Schorgofer, Icarus (2010)
+  if((Jdry - Jsat).le.0) then
+     index_filling = ilay
+     exit
+   endif
+enddo
+
+if(index_filling.eq.1)  depth_filling = mlayer_PEM(1) 
+
+if(index_filling.gt.1) then
+    Jdry_prevlay = (psat_soil(index_filling-1) - pwat_surf)/mlayer_PEM(index_filling-1)
+    Jsat_prevlay = eta(index_filling-1)*dz_psat(index_filling-1)
+   call findroot(Jdry-Jsat,Jdry_prevlay-Jsat_prevlay,mlayer_PEM(index_filling),mlayer_PEM(index_filling-1),depth_filling)
+endif
+
+
+! 2. Compute d_z (eta* d_z(rho)) (last term in Eq. 13 of Schorgofer, Icarus (2010))
+
+
+! 2.0 preliminary: depth to shallowest  ice (discontinuity at interface)
+
+index_firstice = -1
+do ilay = 1,nsoilmx_PEM
+   if (porefill(ilay).le.0.) then
+        index_firstice = ilay  ! first point with ice
+        exit
+    endif
+enddo
+
+! 2.1: now we can computeCompute d_z (eta* d_z(rho))
+
+call deriv1(mlayer_PEM,nsoilmx_PEM,eta,1.,eta(nsoilmx_PEM-1),dz_eta)
+
+if ((index_firstice.gt.0).and.(index_firstice.lt.nsoilmx_PEM-2)) then
+     call deriv1_onesided(index_firstice,mlayer_PEM,nsoilmx_PEM,eta,dz_eta(index_firstice))
+endif
+
+call deriv2_simple(mlayer_PEM,nsoilmx_PEM,psat_soil,psat_surf,psat_bottom,dzz_psat)
+dz_etadz_rho(:) = pvap2rho*(dz_eta(:)*dz_psat(:) + eta(:)*dzz_psat(:))
+
+! 3. Ice table boundary due to geothermal heating
+
+if(index_IS.gt.0) index_geothermal = -1
+if(index_geothermal.lt.0) depth_geothermal = -1.
+if((index_geothermal.gt.0).and.(index_IS.lt.0)) then ! Eq. 21 from Schorfoger, Icarus (2010)
+   index_geothermal = -1
+   do ilay=2,nsoilmx_PEM
+        if (dz_psat(ilay).gt.0.) then  ! first point with reversed flux
+           index_geothermal=ilay
+           call findroot(dz_psat(ilay-1),dz_psat(ilay),mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
+           exit
+        endif
+     enddo
+  else
+     index_geothermal = -1
+  endif
+ if ((index_geothermal.gt.0).and.(index_IS.lt.0)) then !  Eq. 24 from Schorgofer, Icarus (2010)
+     call  colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,index_geothermal-1,nsoilmx_PEM,massfillabove) 
+     index_tmp = -1
+     do ilay=index_geothermal,nsoilmx_PEM
+        if (minval(eta(ilay:nsoilmx_PEM)).le.0.) cycle  ! eta=0 means completely full
+         call colint(porefill(:)/eta(:),mlayer_PEM,nsoilmx_PEM,ilay,nsoilmx_PEM,massfillafter)
+        if (massfillafter<dz_psat(ilay)*pvap2rho*B) then  ! usually executes on i=typeG
+           if (ilay>index_geothermal) then
+!              write(34,*) '# adjustment to geotherm depth by',ilay-index_geothermal
+              call findroot(dz_psat(ilay-1)*pvap2rho*B-massfillabove, &  
+                   dz_psat(ilay)*pvap2rho*B-massfillafter,mlayer_PEM(ilay-1),mlayer_PEM(ilay),depth_geothermal)
+!               if (depth_geothermal.gt.mlayer_PEM(ilay) .or. depth_geothermal.lt.<mlayer_PEM(ilay-1)) write(34,*) 
+!                     '# WARNING: zdepthG interpolation failed',ilay,mlayer_PEM(ilay-1),depth_geothermal,mlayer_PEM(ilay)
+              index_tmp=ilay
+           endif
+           ! otherwise leave depth_geothermal unchanged
+           exit
+        endif
+        massfillabove = massfillafter
+     enddo
+     if (index_tmp.gt.0) index_geothermal = index_tmp
+  end if
+  return
+end subroutine 
+
+
+
+
+
+SUBROUTINE findroot(y1,y2,z1,z2,zr)
+        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        !!!
+        !!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
+        !!! 
+        !!! Author: LL
+        !!!
+        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+        implicit none
+        real,intent(in) :: y1,y2 ! difference between surface water density and at depth  [kg/m^3]
+        real,intent(in) :: z1,z2 ! depth [m]
+        real,intent(out) :: zr ! depth at which we have zero
+        zr = (y1*z2 - y2*z1)/(y1-y2)
+        RETURN
+        end 
+
+SUBROUTINE constriction(porefill,eta)
+
+        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        !!!
+        !!! Purpose: Compute the  constriction of vapor flux by pore ice 
+        !!! 
+        !!! Author: LL
+        !!!
+        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        implicit none
+        real,intent(in) :: porefill ! pore filling fraction
+        real,intent(out) :: eta ! constriction
+
+        !!!
+
+
+        if (porefill.le.0.) eta = 1.
+        if ((porefill.gt.0.) .and.(porefill.lt.1.)) then
+        eta = (1-porefill)**2  ! Hudson et al., JGR, 2009
+        endif
+        if (porefill.le.1.) eta = 0.
+        return
+        end
+
+
+
+subroutine deriv1(z,nz,y,y0,ybot,dzY)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid 
+!!!           
+!!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        implicit none
+
+! first derivative of a function y(z) on irregular grid
+! upper boundary conditions: y(0)=y0
+! lower boundary condition.: yp = ybottom
+  integer, intent(IN) :: nz ! number of layer
+  real, intent(IN) :: z(nz) ! depth layer
+  real, intent(IN) :: y(nz) ! function which needs to be derived
+  real, intent(IN) :: y0,ybot ! boundary conditions
+  real, intent(OUT) :: dzY(nz) ! derivative of y w.r.t depth
+! local
+  integer :: j
+  real :: hm,hp,c1,c2,c3
+
+  hp = z(2)-z(1)
+  c1 = z(1)/(hp*z(2))
+  c2 = 1/z(1) - 1/(z(2)-z(1))
+  c3 = -hp/(z(1)*z(2))
+  dzY(1) = c1*y(2) + c2*y(1) + c3*y0
+  do j=2,nz-1
+     hp = z(j+1)-z(j)
+     hm = z(j)-z(j-1)
+     c1 = +hm/(hp*(z(j+1)-z(j-1)))
+     c2 = 1/hm - 1/hp
+     c3 = -hp/(hm*(z(j+1)-z(j-1)))
+     dzY(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
+  enddo
+  dzY(nz) = (ybot - y(nz-1))/(2.*(z(nz)-z(nz-1)))
+ return
+end subroutine deriv1
+
+
+
+subroutine deriv2_simple(z,nz,y,y0,yNp1,yp2)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid 
+!!!           
+!!! Author: N.S (raw copy/paste from MSIM)
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  ! second derivative y_zz on irregular grid
+  ! boundary conditions: y(0)=y0, y(nz+1)=yNp1
+  implicit none
+  integer, intent(IN) :: nz
+  real, intent(IN) :: z(nz),y(nz),y0,yNp1
+  real, intent(OUT) :: yp2(nz)
+  integer j
+  real hm,hp,c1,c2,c3
+
+  c1 = +2./((z(2)-z(1))*z(2))
+  c2 = -2./((z(2)-z(1))*z(1))
+  c3 = +2./(z(1)*z(2))
+  yp2(1) = c1*y(2) + c2*y(1) + c3*y0
+
+  do j=2,nz-1
+     hp = z(j+1)-z(j)
+     hm = z(j)-z(j-1)
+     c1 = +2./(hp*(z(j+1)-z(j-1)))
+     c2 = -2./(hp*hm)
+     c3 = +2./(hm*(z(j+1)-z(j-1)))
+     yp2(j) = c1*y(j+1) + c2*y(j) + c3*y(j-1)
+  enddo
+  yp2(nz) = (yNp1 - 2*y(nz) + y(nz-1))/(z(nz)-z(nz-1))**2
+  return
+end subroutine deriv2_simple
+
+
+
+subroutine  deriv1_onesided(j,z,nz,y,dy_zj)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose:   First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
+!!!           
+!!! Author: N.S (raw copy/paste from MSIM)
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  implicit none
+  integer, intent(IN) :: nz,j
+  real, intent(IN) :: z(nz),y(nz)
+  real, intent(out) :: dy_zj
+  real h1,h2,c1,c2,c3
+  if (j<1 .or. j>nz-2) then
+     dy_zj = -1.
+  else
+     h1 = z(j+1)-z(j)
+     h2 = z(j+2)-z(j+1)
+     c1 = -(2*h1+h2)/(h1*(h1+h2))
+     c2 =  (h1+h2)/(h1*h2)
+     c3 = -h1/(h2*(h1+h2))
+     dy_zj = c1*y(j) + c2*y(j+1) + c3*y(j+2)
+  endif
+  return
+end subroutine deriv1_onesided
+
+
+subroutine  colint(y,z,nz,i1,i2,integral)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose:  Column integrates y on irregular grid
+!!!           
+!!! Author: N.S (raw copy/paste from MSIM)
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  
+  implicit none
+  integer, intent(IN) :: nz, i1, i2
+  real, intent(IN) :: y(nz), z(nz)
+  real,intent(out) :: integral
+  integer i
+  real dz(nz)
+  
+  dz(1) = (z(2)-0.)/2
+  do i=2,nz-1
+     dz(i) = (z(i+1)-z(i-1))/2.
+  enddo
+  dz(nz) = z(nz)-z(nz-1)
+  integral = sum(y(i1:i2)*dz(i1:i2))
+end subroutine colint
+
+
+
+
+
+
+        end module
Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2901)
+++ /trunk/LMDZ.MARS/README	(revision 2902)
@@ -3904,3 +3904,4 @@
 nslope is now read and written in the startfi.nc
 
-
+== 22/02/2023 == LL
+PEM: First implementation of a dynamic ice table (vs. static ice table at equilibrium) following Schorgofer (2010,Icarus)
Index: /trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F	(revision 2901)
+++ /trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F	(revision 2902)
@@ -170,5 +170,6 @@
       REAL masseq(nlayer),wq(nlayer+1)
       INTEGER ifils,iq2
-
+c LL: flux dependant albedo
+      real :: totflux(ngrid)
 c----------------------------------------------------------------------
 
@@ -176,4 +177,5 @@
 c   --------------
 c
+       totflux(:) = fluxsurf_sw(:,1)+fluxsurf_sw(:,2)
       ! AS: firstcall OK absolute
       IF (firstcall) THEN
@@ -500,4 +502,16 @@
 !     ----------------------------------------
       CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref)
+      DO ig = 1,ngrid
+        if(piceco2(ig).gt.0.) then
+          psolaralb(ig,1) = 0.267+0.00059*totflux(ig)
+          psolaralb(ig,2) = psolaralb(ig,1)
+          if(ig.eq.ngrid) then
+           psolaralb(ig,1) = 0.532+8.72e-4*totflux(ig)
+           psolaralb(ig,2) = psolaralb(ig,1)
+          endif
+        endif
+      ENDDO
+
+
 
 ! set pemisurf() to emissiv when there is bare surface (needed for co2snow)
@@ -739,8 +753,8 @@
 !             write(*,*) "co2condens: South pole: latitude(ngrid)=",
 !     &                                           latitude(ngrid)
-             ztcondsol(ngrid)=
-     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
-     &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
-             pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
+!             ztcondsol(ngrid)=
+!     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*
+!     &                    (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep)))
+!             pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep
            endif
          endif