Changeset 3261 for trunk/LMDZ.MARS/libf

Timestamp:

Mar 12, 2024, 4:37:03 PM (10 months ago)

Author:

llange

Message:

Mars PCM
Reversing -r 3250 (the ""bug"" I fixed was not a bug but my mystake)
LL

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• co2condens_mod.F (modified) (1 diff)
• comsoil_h.F90 (modified) (2 diffs)
• conf_phys.F (modified) (2 diffs)
• dyn1d/init_testphys1d_mod.F90 (modified) (7 diffs)
• dyn1d/testphys1d.F90 (modified) (1 diff)
• physiq_mod.F (modified) (7 diffs)
• soilwater.F90 (modified) (9 diffs)
• vdifc_mod.F (modified) (17 diffs)

trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

@@ -5 +5 @@
       logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2
 !$OMP THREADPRIVATE(scavco2cond)
-      real,    save :: CO2cond_ps = 1.       ! Coefficient to control the surface pressure change
+      real,    save :: CO2cond_ps = 0.       ! Coefficient to control the surface pressure change
 
       CONTAINS

trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90

@@ -31 +31 @@
 ! Subsurface tracers:
 logical, save                                  :: adsorption_soil      ! boolean to call adosrption (or not)
+logical, save                                  :: ads_const_D          ! boolean to have constant diffusion coefficient
+logical, save                                  :: ads_massive_ice      ! boolean to have massive subsurface ice
 real,    save                                  :: choice_ads           ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
 real,    save, allocatable, dimension(:,:,:,:) :: qsoil                ! subsurface tracers (kg/m^3 of regol)
@@ -39 +41 @@
 REAL, parameter :: porosity_reg  = 0.45
 !$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
+!$OMP THREADPRIVATE(ads_const_D,ads_massive_ice)
 
 !=======================================================================

trunk/LMDZ.MARS/libf/phymars/conf_phys.F

@@ -40 +40 @@
      &                           lag_layer 
       use microphys_h, only: mteta
-      use comsoil_h, only: adsorption_soil, choice_ads
+      use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D,
+     & ads_massive_ice
       use nonoro_gwd_mix_mod, only: calljliu_gwimix
 
@@ -1042 +1043 @@
          call getin_p("choice_ads",choice_ads) 
          
+         
+         ads_const_D = .false.
+         call getin_p("ads_const_D",ads_const_D) 
+         
+         ads_massive_ice = .false.
+         call getin_p("ads_massive_ice",ads_massive_ice) 
          poreice_tifeedback = .false.
          call getin_p("poreice_tifeedback",poreice_tifeedback)

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

@@ -22 +22 @@
 use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, &
                                     adsorption_soil, qsoil, igcm_h2o_ice_soil, porosity_reg, &
-                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var
+                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var, ads_massive_ice
 use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
 use dimradmars_mod,           only: tauvis, totcloudfrac, albedo, ini_dimradmars_mod_slope_var, end_dimradmars_mod_slope_var
@@ -116 +116 @@
 ! LL: Subsurface water ice 
 real :: rho_H2O_ice           ! Density (kg/m^3) of water ice
+logical :: ice_fb = .false.
 
 !=======================================================================
@@ -669 +670 @@
 if (.not. therestartfi) qsoil = 0.
 
+call getin("ice_fb_test1D",ice_fb)
 if (.not. therestartfi) then
     ! Initialize depths
@@ -688 +690 @@
             if(adsorption_soil) then
             ! add subsurface ice in qsoil if one runs with adsorption
-                qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*porosity_reg*rho_H2O_ice
-                qsoil(1,2:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
+               if(ads_massive_ice) then
+                   qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*rho_H2O_ice ! assume no porosity in the ice
+                   qsoil(1,2:,igcm_h2o_ice_soil,:) = rho_H2O_ice
+               else
+                   qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*porosity_reg*rho_H2O_ice
+                   qsoil(1,2:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
+               endif
             endif
         else ! searching for the ice/regolith boundary:
             do isoil = 1,nsoil
                 if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
-                    iref = isoil 
+                    iref = isoil +1
                     exit
                 endif
@@ -701 +708 @@
             inertiedat(1,:iref - 1) = inertiedat(1,1)
             ! We compute the transition in layer(iref)
-            inertiedat(1,iref) = sqrt((layer(iref+1) - layer(iref))/(((ice_depth - layer(iref))/inertiedat(1,1)**2) + ((layer(iref+1) - ice_depth)/inertieice**2)))
+            inertiedat(1,iref) =  sqrt((layer(iref) - layer(iref - 1))/(((ice_depth - layer(iref - 1))/inertiedat(1,1)**2) + ((layer(iref) - ice_depth)/inertieice**2)))
             ! Finally, we compute the underlying ice:
             inertiedat(1,iref + 1:) = inertieice
@@ -708 +715 @@
             ! add subsurface ice in qsoil if one runs with adsorption
                 qsoil(1,:iref - 1,igcm_h2o_ice_soil,:) = 0.
-                qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
-                qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref+1))/(layer(iref) - layer(iref+1))*porosity_reg*rho_H2O_ice
+                if(ads_massive_ice) then
+                   qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = rho_H2O_ice
+                   qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref))/(layer(iref-1) - layer(iref))*rho_H2O_ice
+                else
+                   qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
+                   qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref))/(layer(iref-1) - layer(iref))*porosity_reg*rho_H2O_ice                
+                endif
             endif
         endif ! (ice_depth < layer(1))
@@ -715 +727 @@
         inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
     endif ! ice_depth > 0
-    
+
+if(.not.(ice_fb)) then 
+inertiedat(1,:) = inertiedat(1,1)
+end if 
     do isoil = 1,nsoil
         inertiesoil(1,isoil,:) = inertiedat(1,isoil)

trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90

@@ -257 +257 @@
     ! Compute pressure for next time step
     ! -----------------------------------
-    psurf = psurf + dttestphys*dpsurf(1) ! Surface pressure change
+    psurf = psurf! + dttestphys*dpsurf(1) ! Surface pressure change
     plev = ap + psurf*bp
     play = aps + psurf*bps

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -82 +82 @@
       USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
       USE calldrag_noro_mod, ONLY: calldrag_noro
-      USE vdifc_mod, ONLY: vdifc
+      USE vdifc_mod, ONLY: vdifc, compute_Tice
       use param_v4_h, only: nreact,n_avog,
      &                      fill_data_thermos, allocate_param_thermos
@@ -124 +124 @@
       use lmdz_atke_turbulence_ini, only : atke_ini
       use waterice_tifeedback_mod, only : waterice_tifeedback
+      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
       IMPLICIT NONE
 c=======================================================================
@@ -542 +543 @@
 ! Variable for ice table
       REAL :: rhowater_surf(ngrid,nslope)            ! Water density at the surface [kg/m^3]
+      REAL :: rhowater_surf_schorgo(ngrid,nslope)    ! Water density at the surface, schorghofer way [kg/m^3]
       REAL :: rhowater_surf_sat(ngrid,nslope)        ! Water density at the surface at saturation [kg/m^3]
       REAL :: rhowater_soil(ngrid,nsoilmx,nslope)    ! Water density in soil layers [kg/m^3]
+      REAL :: rhowater_ice(ngrid,nslope)             ! Water density at the subsurface ice depth[kg/m^3]
       REAL,PARAMETER :: alpha_clap_h2o = 28.9074     ! Coeff for Clapeyron law [/]
       REAL,PARAMETER :: beta_clap_h2o = -6143.7      ! Coeff for Clapeyron law [K]
@@ -551 +554 @@
       REAL :: ztmp1,ztmp2                            ! intermediate variables to compute the mean molar mass of the layer
       REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
+      REAL :: Tice
 ! Variable for the computation of the TKE with parameterization from ATKE working group
       REAL :: viscom                              ! kinematic molecular viscosity for momentum
@@ -3939 +3943 @@
      &         * mmol(igcm_h2o_vap)/(mugaz*r)
              endif
+ ! schorghofer way
+             rhowater_surf_schorgo(ig,islope) = min(pvap_surf(ig),
+     &         exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o)
+     &         / tsurf(ig,islope))/tsurf(ig,islope)
+     &         * mmol(igcm_h2o_vap)/(mugaz*r)
              DO isoil = 1,nsoilmx
              rhowater_soil(ig,isoil,islope) =
@@ -3945 +3954 @@
      &         * mmol(igcm_h2o_vap)/(mugaz*r)
              ENDDO
+              if(h2o_ice_depth(ig,islope).gt.0) then
+              call compute_Tice(nsoilmx,tsoil(ig,:,islope),
+     &                     tsurf(ig,islope),h2o_ice_depth(ig,islope),
+     &                     Tice)
+              rhowater_ice(ig,islope) =  
+     &         exp(beta_clap_h2o/Tice+alpha_clap_h2o)
+     &         / Tice
+     &         * mmol(igcm_h2o_vap)/(mugaz*r)
+              else
+              rhowater_ice(ig,islope) = 0.
+              endif
+     
           ENDDO
         ENDDO
@@ -3954 +3975 @@
      &     "rhowater_surface",'kg.m-3',
      &     rhowater_surf(:,iflat))
+     
+           CALL write_output("waterdensity_surface_schorgho",
+     &     "rhowater_surf_schorgo",'kg.m-3',
+     &     rhowater_surf_schorgo(:,iflat))
+     
+                CALL write_output("waterdensity_ssice",
+     &     "rhowater_ice",'kg.m-3',
+     &     rhowater_ice(:,iflat))
       DO islope = 1,nslope
         write(str2(1:2),'(i2.2)') islope

trunk/LMDZ.MARS/libf/phymars/soilwater.F90

@@ -1 +1 @@
 subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep,  & 
       exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf,  & 
-      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice) 
-
-
-      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg
+      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice,var_flux_soil) 
+
+
+      use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg,ads_const_D
       use comcstfi_h
       use tracer_mod
@@ -10 +10 @@
       use geometry_mod, only: cell_area, latitude_deg
       use write_output_mod, only: write_output
+      use paleoclimate_mod, only: d_coef
 implicit none 
 
@@ -183 +184 @@
 logical, allocatable, save :: close_top(:, :)         ! closing diffusion at the top of a layer if ice rises over saturation
 logical, allocatable, save :: close_bottom(:, :)      ! closing diffusion at the bottom of a layer if ice risies over saturation
-logical, parameter :: print_closure_warnings = .true.
+logical, parameter :: print_closure_warnings = .false.
 
 ! Coefficients for the Diffusion calculations
@@ -269 +270 @@
 real*8 :: diff(ngrid, nsoil)              ! difference between total water content and ice + vapor content 
                                           ! (only used for output, inconstistent: should just be adswater)
-real :: var_flux_soil(ngrid, nsoil)
+real,intent(out) :: var_flux_soil(ngrid, nsoil)
 
 ! Loop variables and counters
@@ -568 +569 @@
                   ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 -  only exact for normal diffusion)
                   D_mid(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice(ig, ik)) ** 2.D0 * 1.D0 / (1.D0 / Dm(ig, ik) + 1.D0 / Dk(ig, ik))
-
+                  if(ads_const_D) D_mid(ig, ik) = d_coef(ig,1)
                   ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991)
                   !                             D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik)  &
@@ -593 +594 @@
                   
                   saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
-                  
+                 
                   D_inter(ig, ik) = D0(ig, ik) * (1.D0 - saturation_water_ice_inter(ig, ik)) ** 2.D0 * 1.D0 / (1.D0 / Dm_inter(ig, ik) + 1.D0 / Dk_inter(ig, ik))
+                  
+                  if(ads_const_D) D_inter(ig, ik) = d_coef(ig,1)
+                  
                   
 !                  D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) &  ! old implementation with direct interpolation
@@ -1567 +1571 @@
 !       call write_output("dztot1","change in dztot", "unknown",real(dztot1))
 
-        call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil)                                
       
 !       call write_output("A_coef","A coefficient of dztot1", "unknown",real(B))                                
@@ -1595 +1598 @@
        call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater))                   
       
-!       call write_output("subsurfice","subsurface ice", "kg/m^3",real(ice))                   
+        call write_output("subsurfice_mass_py","subsurface ice", "kg/m^3",real(ice))                   
 
         call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego)                   
@@ -1603 +1606 @@
 !       call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out))                         
 
-!       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter)                         
+       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_mid)                         
             
 !endif

trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

@@ -35 +35 @@
       use microphys_h, only: To
       use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
-      use comsoil_h, only: layer, mlayer,adsorption_soil
+      use comsoil_h, only: layer, mlayer,adsorption_soil,
+     &                     igcm_h2o_ice_soil
       use vdif_cd_mod, only: vdif_cd
       use lmdz_call_atke, only: call_atke
@@ -234 +235 @@
       LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
       REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores
+
+      REAL :: zdqsoildif_tot(ngrid,nsoil,nqsoil,nslope)
+      REAL :: zqsoil(ngrid,nsoil,nqsoil,nslope)
 !! Water buyoncy
       LOGICAL :: virtual
@@ -246 +250 @@
       REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep
       REAL zdqsdif_ssi_tot(ngrid,nslope)       ! Total flux of the interactions with SSI kg/m^2/s^-1)
+      REAL zdqsdif_regolith_frost_tot(ngrid,nslope)
+      REAL zdqsdif_regolith_atm_tot(ngrid,nslope)
+      real flux_soil_py(ngrid,nsoil,nslope)
+      real flux_soil_py_tot(ngrid,nsoil,nslope)
+      
+      real flux_schorgho(ngrid,nslope)
+      real flux_schorgho_vfrost(ngrid,nslope)
+      
+      real flux_schorgho_tot(ngrid,nslope)
+      real flux_schorgho_vfrost_tot(ngrid,nslope)      
       
       REAL :: tol_frost = 1e-4 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect
@@ -306 +320 @@
       dwatercap_dif(1:ngrid,1:nslope)=0
       zdqsdif_regolith(1:ngrid,1:nslope)=0
+      zdqsdif_regolith_frost_tot(1:ngrid,1:nslope)=0
+      zdqsdif_regolith_atm_tot(1:ngrid,1:nslope)=0
       zq1temp_regolith(1:ngrid)=0 
       zdqsdif_tot(1:ngrid)=0
@@ -315 +331 @@
       zdqsdif_ssi_atm_tot(:,:) = 0.
       zdqsdif_ssi_frost_tot(:,:) = 0.
-      
-           
+      zdqsoildif_tot(:,:,:,:) = 0.
+      zqsoil(:,:,:,:) = 0. 
+      flux_soil_py(:,:,:) = 0.
+      flux_soil_py_tot(:,:,:) = 0.     
+      flux_schorgho(:,:) = 0.
+      flux_schorgho_vfrost(:,:) = 0.
+      flux_schorgho_tot(:,:) = 0.
+      flux_schorgho_vfrost_tot(:,:)  = 0.
+      
+      
 c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
 c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
@@ -1004 +1028 @@
               watercap_tmp(ig) = watercap(ig,islope)/
      &                       cos(pi*def_slope_mean(islope)/180.)
+
+              zqsoil(ig,:,:,islope) = qsoil(ig,:,:,islope)/
+     &             cos(pi*def_slope_mean(islope)/180.)
+
            ENDDO ! ig=1,ngrid
            
@@ -1015 +1043 @@
            saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)    
            DO ig=1,ngrid
-c           nsubtimestep(ig)=1 !for debug
+!           nsubtimestep(ig)=1 !for debug
            subtimestep = ptimestep/nsubtimestep(ig)
                           call write_output('subtimestep',
@@ -1089 +1117 @@
      &                     za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:), 
      &                     zdqsdif_surf(ig), zqsurf(ig), 
-     &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig), 
+     &                     zqsoil(ig,:,:,islope), pplev(ig,1), rho(ig), 
      &                     writeoutput,zdqsdif_regolith(ig,islope), 
      &                     zq1temp_regolith(ig),
-     &                     pore_icefraction(ig,:,islope))
-
+     &                     pore_icefraction(ig,:,islope),
+     &                     flux_soil_py(ig,:,islope))
+                flux_soil_py_tot(ig,:,islope) = 
+     &           flux_soil_py_tot(ig,:,islope) +
+     &           flux_soil_py(ig,:,islope)*subtimestep
 
                  if(.not.watercaptag(ig)) then
@@ -1100 +1131 @@
                          zdqsdif_tot(ig)=
      &                        -zqsurf(ig)/subtimestep
+                         zdqsdif_regolith_atm_tot(ig,islope) =
+     &                         zdqsdif_regolith_atm_tot(ig,islope) 
+     &                   + zdqsdif_regolith(ig,islope)*subtimestep
                      else
                          zdqsdif_tot(ig) = zdqsdif_surf(ig) +
      &                           zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface)
+                          zdqsdif_regolith_frost_tot(ig,islope) = 
+     &                     zdqsdif_regolith_frost_tot(ig,islope) + 
+     &                     zdqsdif_regolith(ig,islope)*subtimestep
                      endif  ! of "if exchange = true"
                  endif  ! of "if not.watercaptag"
@@ -1146 +1183 @@
                          qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
      &                                       *qsat_ssi(ig,islope)
+
+                  ! And now we solve correctly the system      
+                        z1(ig)=1./(za(ig,1)+zb(ig,1)+
+     &                         zb(ig,2)*(1.-zd(ig,2)))
+                        zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
+     &                           zb(ig,2)*zc(ig,2) + 
+     &                           (-zdqsdif_tot(ig)) *subtimestep)
+     &                           * z1(ig)
+                        zd(ig,1)=zb(ig,1)*z1(ig)
+                        zq1temp(ig)=zc(ig,1)+ zd(ig,1)
+     &                              *qsat_ssi(ig,islope)
+            
                   ! Flux from the subsurface      
                         if(old_wsublimation_scheme) then
@@ -1158 +1207 @@
      &                     *(zq1temp(ig)-qsat_ssi(ig,islope))                         
                         endif
-                  ! And now we solve correctly the system      
-                        z1(ig)=1./(za(ig,1)+zb(ig,1)+
-     &                         zb(ig,2)*(1.-zd(ig,2)))
-                        zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
-     &                           zb(ig,2)*zc(ig,2) + 
-     &                           (-zdqsdif_tot(ig)) *subtimestep)
-     &                           * z1(ig)
-                        zd(ig,1)=zb(ig,1)*z1(ig)
-                        zq1temp(ig)=zc(ig,1)+ zd(ig,1)
-     &                              *qsat_ssi(ig,islope)
-            
-
 
                     else 
@@ -1179 +1216 @@
      &                      zb(ig,2)*zc(ig,2) +
      &                      (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
-                            zq1temp(ig)=zc(ig,1)
+                         zq1temp(ig)=zc(ig,1)
                     endif ! Of subsurface <-> atmosphere exchange
                  endif ! sublimating more than available frost & surface - frost exchange 
@@ -1214 +1251 @@
      &                                 *(qsat(ig)-qsat_ssi(ig,islope))
 
-! Line to comment for now since we don't have a mass of subsurface frost in our computation (otherwise, we would not conserve the H2O mass in the system)               
+! Line to comment for now since we don't have a mass of subsurface frost in our computation (otherwise, we would not conserve the H2O mass in the system)              
+
+                  if ((-zdqsdif_tot(ig)+zdqsdif_ssi_frost(ig,islope))
+     &            *subtimestep.ge.(zqsurf(ig))) then
+                 
+                     zdqsdif_ssi_frost(ig,islope) = 
+     &                         zqsurf(ig)/subtimestep
+     &                         +   zdqsdif_tot(ig)
+                  endif
+                  
                     zdqsdif_tot(ig) = zdqsdif_tot(ig) - 
-     &                        zdqsdif_ssi_frost(ig,islope) 
+     &                        zdqsdif_ssi_frost(ig,islope)
+
                  endif ! watercaptag or frost at the surface
              endif ! interaction frost <-> subsurface ice
              
 
+ccc Special case : schorgho study
+
+       if (h2o_ice_depth(ig,islope).gt.0) then
+                    call watersat(1,Tice(ig,islope),pplev(ig,1),
+     &                            qsat_ssi(ig,islope))
+                  qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
+     &                                       *qsat_ssi(ig,islope)
+           flux_schorgho(ig,islope) = d_coef(ig,islope)
+     &                               /h2o_ice_depth(ig,islope)
+     &                              *rho(ig)*dryness(ig)
+     &                              *(min(zq1temp(ig),qsat(ig))
+     &                               -qsat_ssi(ig,islope))
+     
+            if((zqsurf(ig).gt.tol_frost)) then
+            flux_schorgho_vfrost(ig,islope) = d_coef(ig,islope)
+     &                               /h2o_ice_depth(ig,islope)
+     &                              *rho(ig)*dryness(ig)
+     &                              *(qsat(ig)
+     &                               -qsat_ssi(ig,islope))
+             else
+            flux_schorgho_vfrost(ig,islope) = d_coef(ig,islope)
+     &                               /h2o_ice_depth(ig,islope)
+     &                              *rho(ig)*dryness(ig)
+     &                              *(zq1temp(ig)
+     &                               -qsat_ssi(ig,islope))
+                 
+             endif
+           
+        endif
+
+ccc
 c             Starting upward calculations for water :
              zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
@@ -1246 +1324 @@
      &                              zdqsdif_ssi_atm_tot(ig,islope)
      &                            + zdqsdif_ssi_atm(ig,islope)
+     &                            * subtimestep
               zdqsdif_ssi_frost_tot(ig,islope) = 
      &                                zdqsdif_ssi_frost_tot(ig,islope) 
      &                              + zdqsdif_ssi_frost(ig,islope)
+     &                              * subtimestep
+     
+              flux_schorgho_vfrost_tot(ig,islope) = 
+     &           flux_schorgho_vfrost_tot(ig,islope) + 
+     &           flux_schorgho_vfrost(ig,islope)* subtimestep
+
+              flux_schorgho_tot(ig,islope) = 
+     &           flux_schorgho_tot(ig,islope) + 
+     &           flux_schorgho(ig,islope)* subtimestep     
+     
+     
 c             Monitoring instantaneous latent heat flux in W.m-2 :
               zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
@@ -1276 +1366 @@
      &                       cos(pi*def_slope_mean(islope)/180.))
      &                       /ptimestep
-     
+    
+
+            zdqsoildif_tot(ig,:,:,islope) =  (zqsoil(ig,:,:,islope) - 
+     &                      qsoil(ig,:,:,islope)/
+     &                       cos(pi*def_slope_mean(islope)/180.))
+     &                       /ptimestep
+            qsoil(ig,:,:,islope) = zqsoil(ig,:,:,islope)
             zdqsdif_ssi_tot(ig,islope) = 
      &                        zdqsdif_ssi_atm_tot(ig,islope) 
@@ -1413 +1509 @@
      &                zdqsdif_ssi_tot(:,iflat))
 
-
+        call write_output('zdqsdif_regolith_frost_tot',
+     &                '','',zdqsdif_regolith_frost_tot(:,iflat))
+
+        call write_output('zdqsdif_regolith_atm_tot',
+     &                '','',zdqsdif_regolith_atm_tot(:,iflat))
+
+
+        call write_output('zdqsoildif_tot','','',
+     &           zdqsoildif_tot(:,:,igcm_h2o_ice_soil,iflat))
+     
+             call write_output('zdqsoildif_tot','','',
+     &           zdqsoildif_tot(:,:,igcm_h2o_ice_soil,iflat))
+    
+    
+         call write_output('flux_schorgho_vfrost_tot',
+     &     '','kg.m-2.s-1',
+     &                flux_schorgho_vfrost_tot(:,iflat))
+     
+         call write_output('flux_schorgho_tot',
+     &     '','kg.m-2.s-1',
+     &                flux_schorgho_tot(:,iflat)) 
 C       Diagnostic output for HDO
 !        if (hdo) then
@@ -1423 +1539 @@
 !     &                       ' ',h2oflux(:))
 !        endif
+
+         call write_output("flux_soillayer","",
+     &             "",flux_soil_py_tot(:,:,iflat))
+
 
 c-----------------------------------------------------------------------