Changeset 3261 for trunk


Ignore:
Timestamp:
Mar 12, 2024, 4:37:03 PM (9 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
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3253 r3261  
    45454545Following -r 3098; cleaning of vdfic for the management of subsurface water ice.
    45464546Fixing some errors (wrong interpolation to compute the water ice temperature, wrong boundary conditions to compute qvap(1))
     4547
     4548== 12/03/2024 == LL
     4549Reversing -r 3250 (the ""bug"" I fixed was not a bug but my mystake)
     4550

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

    r3207 r3261  
    55      logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2
    66!$OMP THREADPRIVATE(scavco2cond)
    7       real,    save :: CO2cond_ps = 1.       ! Coefficient to control the surface pressure change
     7      real,    save :: CO2cond_ps = 0.       ! Coefficient to control the surface pressure change
    88
    99      CONTAINS

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

    r3230 r3261  
    3131! Subsurface tracers:
    3232logical, save                                  :: adsorption_soil      ! boolean to call adosrption (or not)
     33logical, save                                  :: ads_const_D          ! boolean to have constant diffusion coefficient
     34logical, save                                  :: ads_massive_ice      ! boolean to have massive subsurface ice
    3335real,    save                                  :: choice_ads           ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
    3436real,    save, allocatable, dimension(:,:,:,:) :: qsoil                ! subsurface tracers (kg/m^3 of regol)
    … …  
    3941REAL, parameter :: porosity_reg  = 0.45
    4042!$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
     43!$OMP THREADPRIVATE(ads_const_D,ads_massive_ice)
    4144
    4245!=======================================================================

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

    r3230 r3261  
    4040     &                           lag_layer
    4141      use microphys_h, only: mteta
    42       use comsoil_h, only: adsorption_soil, choice_ads
     42      use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D,
     43     & ads_massive_ice
    4344      use nonoro_gwd_mix_mod, only: calljliu_gwimix
    4445
    … …  
    10421043         call getin_p("choice_ads",choice_ads)
    10431044         
     1045         
     1046         ads_const_D = .false.
     1047         call getin_p("ads_const_D",ads_const_D)
     1048         
     1049         ads_massive_ice = .false.
     1050         call getin_p("ads_massive_ice",ads_massive_ice)
    10441051         poreice_tifeedback = .false.
    10451052         call getin_p("poreice_tifeedback",poreice_tifeedback)

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

    r3253 r3261  
    2222use comsoil_h,                only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, &
    2323                                    adsorption_soil, qsoil, igcm_h2o_ice_soil, porosity_reg, &
    24                                     ini_comsoil_h_slope_var, end_comsoil_h_slope_var
     24                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var, ads_massive_ice
    2525use comvert_mod,              only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight
    2626use dimradmars_mod,           only: tauvis, totcloudfrac, albedo, ini_dimradmars_mod_slope_var, end_dimradmars_mod_slope_var
    … …  
    116116! LL: Subsurface water ice
    117117real :: rho_H2O_ice           ! Density (kg/m^3) of water ice
     118logical :: ice_fb = .false.
    118119
    119120!=======================================================================
    … …  
    669670if (.not. therestartfi) qsoil = 0.
    670671
     672call getin("ice_fb_test1D",ice_fb)
    671673if (.not. therestartfi) then
    672674    ! Initialize depths
    … …  
    688690            if(adsorption_soil) then
    689691            ! add subsurface ice in qsoil if one runs with adsorption
    690                 qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*porosity_reg*rho_H2O_ice
    691                 qsoil(1,2:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
     692               if(ads_massive_ice) then
     693                   qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*rho_H2O_ice ! assume no porosity in the ice
     694                   qsoil(1,2:,igcm_h2o_ice_soil,:) = rho_H2O_ice
     695               else
     696                   qsoil(1,1,igcm_h2o_ice_soil,:) = (ice_depth-layer(2))/(layer(1) - layer(2))*porosity_reg*rho_H2O_ice
     697                   qsoil(1,2:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
     698               endif
    692699            endif
    693700        else ! searching for the ice/regolith boundary:
    694701            do isoil = 1,nsoil
    695702                if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
    696                     iref = isoil
     703                    iref = isoil +1
    697704                    exit
    698705                endif
    … …  
    701708            inertiedat(1,:iref - 1) = inertiedat(1,1)
    702709            ! We compute the transition in layer(iref)
    703             inertiedat(1,iref) = sqrt((layer(iref+1) - layer(iref))/(((ice_depth - layer(iref))/inertiedat(1,1)**2) + ((layer(iref+1) - ice_depth)/inertieice**2)))
     710            inertiedat(1,iref) =  sqrt((layer(iref) - layer(iref - 1))/(((ice_depth - layer(iref - 1))/inertiedat(1,1)**2) + ((layer(iref) - ice_depth)/inertieice**2)))
    704711            ! Finally, we compute the underlying ice:
    705712            inertiedat(1,iref + 1:) = inertieice
    … …  
    708715            ! add subsurface ice in qsoil if one runs with adsorption
    709716                qsoil(1,:iref - 1,igcm_h2o_ice_soil,:) = 0.
    710                 qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
    711                 qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref+1))/(layer(iref) - layer(iref+1))*porosity_reg*rho_H2O_ice
     717                if(ads_massive_ice) then
     718                   qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = rho_H2O_ice
     719                   qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref))/(layer(iref-1) - layer(iref))*rho_H2O_ice
     720                else
     721                   qsoil(1,iref+1:,igcm_h2o_ice_soil,:) = porosity_reg*rho_H2O_ice
     722                   qsoil(1,iref,igcm_h2o_ice_soil,:) = (ice_depth-layer(iref))/(layer(iref-1) - layer(iref))*porosity_reg*rho_H2O_ice               
     723                endif
    712724            endif
    713725        endif ! (ice_depth < layer(1))
    … …  
    715727        inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
    716728    endif ! ice_depth > 0
    717    
     729
     730if(.not.(ice_fb)) then
     731inertiedat(1,:) = inertiedat(1,1)
     732end if
    718733    do isoil = 1,nsoil
    719734        inertiesoil(1,isoil,:) = inertiedat(1,isoil)

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

    r3231 r3261  
    257257    ! Compute pressure for next time step
    258258    ! -----------------------------------
    259     psurf = psurf + dttestphys*dpsurf(1) ! Surface pressure change
     259    psurf = psurf! + dttestphys*dpsurf(1) ! Surface pressure change
    260260    plev = ap + psurf*bp
    261261    play = aps + psurf*bps

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

    r3230 r3261  
    8282      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
    8383      USE calldrag_noro_mod, ONLY: calldrag_noro
    84       USE vdifc_mod, ONLY: vdifc
     84      USE vdifc_mod, ONLY: vdifc, compute_Tice
    8585      use param_v4_h, only: nreact,n_avog,
    8686     &                      fill_data_thermos, allocate_param_thermos
    … …  
    124124      use lmdz_atke_turbulence_ini, only : atke_ini
    125125      use waterice_tifeedback_mod, only : waterice_tifeedback
     126      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
    126127      IMPLICIT NONE
    127128c=======================================================================
    … …  
    542543! Variable for ice table
    543544      REAL :: rhowater_surf(ngrid,nslope)            ! Water density at the surface [kg/m^3]
     545      REAL :: rhowater_surf_schorgo(ngrid,nslope)    ! Water density at the surface, schorghofer way [kg/m^3]
    544546      REAL :: rhowater_surf_sat(ngrid,nslope)        ! Water density at the surface at saturation [kg/m^3]
    545547      REAL :: rhowater_soil(ngrid,nsoilmx,nslope)    ! Water density in soil layers [kg/m^3]
     548      REAL :: rhowater_ice(ngrid,nslope)             ! Water density at the subsurface ice depth[kg/m^3]
    546549      REAL,PARAMETER :: alpha_clap_h2o = 28.9074     ! Coeff for Clapeyron law [/]
    547550      REAL,PARAMETER :: beta_clap_h2o = -6143.7      ! Coeff for Clapeyron law [K]
    … …  
    551554      REAL :: ztmp1,ztmp2                            ! intermediate variables to compute the mean molar mass of the layer
    552555      REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
     556      REAL :: Tice
    553557! Variable for the computation of the TKE with parameterization from ATKE working group
    554558      REAL :: viscom                              ! kinematic molecular viscosity for momentum
    … …  
    39393943     &         * mmol(igcm_h2o_vap)/(mugaz*r)
    39403944             endif
     3945 ! schorghofer way
     3946             rhowater_surf_schorgo(ig,islope) = min(pvap_surf(ig),
     3947     &         exp(beta_clap_h2o/tsurf(ig,islope)+alpha_clap_h2o)
     3948     &         / tsurf(ig,islope))/tsurf(ig,islope)
     3949     &         * mmol(igcm_h2o_vap)/(mugaz*r)
    39413950             DO isoil = 1,nsoilmx
    39423951             rhowater_soil(ig,isoil,islope) =
    … …  
    39453954     &         * mmol(igcm_h2o_vap)/(mugaz*r)
    39463955             ENDDO
     3956              if(h2o_ice_depth(ig,islope).gt.0) then
     3957              call compute_Tice(nsoilmx,tsoil(ig,:,islope),
     3958     &                     tsurf(ig,islope),h2o_ice_depth(ig,islope),
     3959     &                     Tice)
     3960              rhowater_ice(ig,islope) = 
     3961     &         exp(beta_clap_h2o/Tice+alpha_clap_h2o)
     3962     &         / Tice
     3963     &         * mmol(igcm_h2o_vap)/(mugaz*r)
     3964              else
     3965              rhowater_ice(ig,islope) = 0.
     3966              endif
     3967     
    39473968          ENDDO
    39483969        ENDDO
    … …  
    39543975     &     "rhowater_surface",'kg.m-3',
    39553976     &     rhowater_surf(:,iflat))
     3977     
     3978           CALL write_output("waterdensity_surface_schorgho",
     3979     &     "rhowater_surf_schorgo",'kg.m-3',
     3980     &     rhowater_surf_schorgo(:,iflat))
     3981     
     3982                CALL write_output("waterdensity_ssice",
     3983     &     "rhowater_ice",'kg.m-3',
     3984     &     rhowater_ice(:,iflat))
    39563985      DO islope = 1,nslope
    39573986        write(str2(1:2),'(i2.2)') islope

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

    r3248 r3261  
    11subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep,  &
    22      exchange, qsat_surf, pq, pa, pb, pc, pd, pdqsdifpot, pqsurf,  &
    3       pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice)
    4 
    5 
    6       use comsoil_h, only: igcm_h2o_vap_soil, igcm_h2o_ice_soil, igcm_h2o_vap_ads, layer, mlayer, choice_ads, porosity_reg
     3      pqsoil, pplev, rhoatmo, writeoutput, zdqsdifrego, zq1temp2, saturation_water_ice,var_flux_soil)
     4
     5
     6      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
    77      use comcstfi_h
    88      use tracer_mod
    … …  
    1010      use geometry_mod, only: cell_area, latitude_deg
    1111      use write_output_mod, only: write_output
     12      use paleoclimate_mod, only: d_coef
    1213implicit none
    1314
    … …  
    183184logical, allocatable, save :: close_top(:, :)         ! closing diffusion at the top of a layer if ice rises over saturation
    184185logical, allocatable, save :: close_bottom(:, :)      ! closing diffusion at the bottom of a layer if ice risies over saturation
    185 logical, parameter :: print_closure_warnings = .true.
     186logical, parameter :: print_closure_warnings = .false.
    186187
    187188! Coefficients for the Diffusion calculations
    … …  
    269270real*8 :: diff(ngrid, nsoil)              ! difference between total water content and ice + vapor content
    270271                                          ! (only used for output, inconstistent: should just be adswater)
    271 real :: var_flux_soil(ngrid, nsoil)
     272real,intent(out) :: var_flux_soil(ngrid, nsoil)
    272273
    273274! Loop variables and counters
    … …  
    568569                  ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 -  only exact for normal diffusion)
    569570                  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))
    570 
     571                  if(ads_const_D) D_mid(ig, ik) = d_coef(ig,1)
    571572                  ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991)
    572573                  !                             D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik)  &
    … …  
    593594                 
    594595                  saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
    595                   
     596                 
    596597                  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))
     598                 
     599                  if(ads_const_D) D_inter(ig, ik) = d_coef(ig,1)
     600                 
    597601                 
    598602!                  D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) &  ! old implementation with direct interpolation
    … …  
    15671571!       call write_output("dztot1","change in dztot", "unknown",real(dztot1))
    15681572
    1569         call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil)                               
    15701573     
    15711574!       call write_output("A_coef","A coefficient of dztot1", "unknown",real(B))                               
    … …  
    15951598       call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater))                   
    15961599     
    1597 !       call write_output("subsurfice","subsurface ice", "kg/m^3",real(ice))                   
     1600        call write_output("subsurfice_mass_py","subsurface ice", "kg/m^3",real(ice))                   
    15981601
    15991602        call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego)                   
    … …  
    16031606!       call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out))                         
    16041607
    1605 !       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter)                         
     1608       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_mid)                         
    16061609           
    16071610!endif

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

    r3253 r3261  
    3535      use microphys_h, only: To
    3636      use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer
    37       use comsoil_h, only: layer, mlayer,adsorption_soil
     37      use comsoil_h, only: layer, mlayer,adsorption_soil,
     38     &                     igcm_h2o_ice_soil
    3839      use vdif_cd_mod, only: vdif_cd
    3940      use lmdz_call_atke, only: call_atke
    … …  
    234235      LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
    235236      REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores
     237
     238      REAL :: zdqsoildif_tot(ngrid,nsoil,nqsoil,nslope)
     239      REAL :: zqsoil(ngrid,nsoil,nqsoil,nslope)
    236240!! Water buyoncy
    237241      LOGICAL :: virtual
    … …  
    246250      REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep
    247251      REAL zdqsdif_ssi_tot(ngrid,nslope)       ! Total flux of the interactions with SSI kg/m^2/s^-1)
     252      REAL zdqsdif_regolith_frost_tot(ngrid,nslope)
     253      REAL zdqsdif_regolith_atm_tot(ngrid,nslope)
     254      real flux_soil_py(ngrid,nsoil,nslope)
     255      real flux_soil_py_tot(ngrid,nsoil,nslope)
     256     
     257      real flux_schorgho(ngrid,nslope)
     258      real flux_schorgho_vfrost(ngrid,nslope)
     259     
     260      real flux_schorgho_tot(ngrid,nslope)
     261      real flux_schorgho_vfrost_tot(ngrid,nslope)     
    248262     
    249263      REAL :: tol_frost = 1e-4 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect
    … …  
    306320      dwatercap_dif(1:ngrid,1:nslope)=0
    307321      zdqsdif_regolith(1:ngrid,1:nslope)=0
     322      zdqsdif_regolith_frost_tot(1:ngrid,1:nslope)=0
     323      zdqsdif_regolith_atm_tot(1:ngrid,1:nslope)=0
    308324      zq1temp_regolith(1:ngrid)=0
    309325      zdqsdif_tot(1:ngrid)=0
    … …  
    315331      zdqsdif_ssi_atm_tot(:,:) = 0.
    316332      zdqsdif_ssi_frost_tot(:,:) = 0.
    317      
    318            
     333      zdqsoildif_tot(:,:,:,:) = 0.
     334      zqsoil(:,:,:,:) = 0.
     335      flux_soil_py(:,:,:) = 0.
     336      flux_soil_py_tot(:,:,:) = 0.     
     337      flux_schorgho(:,:) = 0.
     338      flux_schorgho_vfrost(:,:) = 0.
     339      flux_schorgho_tot(:,:) = 0.
     340      flux_schorgho_vfrost_tot(:,:)  = 0.
     341     
     342     
    319343c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
    320344c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
    … …  
    10041028              watercap_tmp(ig) = watercap(ig,islope)/
    10051029     &                       cos(pi*def_slope_mean(islope)/180.)
     1030
     1031              zqsoil(ig,:,:,islope) = qsoil(ig,:,:,islope)/
     1032     &             cos(pi*def_slope_mean(islope)/180.)
     1033
    10061034           ENDDO ! ig=1,ngrid
    10071035           
    … …  
    10151043           saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)   
    10161044           DO ig=1,ngrid
    1017 c           nsubtimestep(ig)=1 !for debug
     1045!           nsubtimestep(ig)=1 !for debug
    10181046           subtimestep = ptimestep/nsubtimestep(ig)
    10191047                          call write_output('subtimestep',
    … …  
    10891117     &                     za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:),
    10901118     &                     zdqsdif_surf(ig), zqsurf(ig),
    1091      &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
     1119     &                     zqsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
    10921120     &                     writeoutput,zdqsdif_regolith(ig,islope),
    10931121     &                     zq1temp_regolith(ig),
    1094      &                     pore_icefraction(ig,:,islope))
    1095 
     1122     &                     pore_icefraction(ig,:,islope),
     1123     &                     flux_soil_py(ig,:,islope))
     1124                flux_soil_py_tot(ig,:,islope) =
     1125     &           flux_soil_py_tot(ig,:,islope) +
     1126     &           flux_soil_py(ig,:,islope)*subtimestep
    10961127
    10971128                 if(.not.watercaptag(ig)) then
    … …  
    11001131                         zdqsdif_tot(ig)=
    11011132     &                        -zqsurf(ig)/subtimestep
     1133                         zdqsdif_regolith_atm_tot(ig,islope) =
     1134     &                         zdqsdif_regolith_atm_tot(ig,islope)
     1135     &                   + zdqsdif_regolith(ig,islope)*subtimestep
    11021136                     else
    11031137                         zdqsdif_tot(ig) = zdqsdif_surf(ig) +
    11041138     &                           zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface)
     1139                          zdqsdif_regolith_frost_tot(ig,islope) =
     1140     &                     zdqsdif_regolith_frost_tot(ig,islope) +
     1141     &                     zdqsdif_regolith(ig,islope)*subtimestep
    11051142                     endif  ! of "if exchange = true"
    11061143                 endif  ! of "if not.watercaptag"
    … …  
    11461183                         qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
    11471184     &                                       *qsat_ssi(ig,islope)
     1185
     1186                  ! And now we solve correctly the system     
     1187                        z1(ig)=1./(za(ig,1)+zb(ig,1)+
     1188     &                         zb(ig,2)*(1.-zd(ig,2)))
     1189                        zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
     1190     &                           zb(ig,2)*zc(ig,2) +
     1191     &                           (-zdqsdif_tot(ig)) *subtimestep)
     1192     &                           * z1(ig)
     1193                        zd(ig,1)=zb(ig,1)*z1(ig)
     1194                        zq1temp(ig)=zc(ig,1)+ zd(ig,1)
     1195     &                              *qsat_ssi(ig,islope)
     1196           
    11481197                  ! Flux from the subsurface     
    11491198                        if(old_wsublimation_scheme) then
    … …  
    11581207     &                     *(zq1temp(ig)-qsat_ssi(ig,islope))                         
    11591208                        endif
    1160                   ! And now we solve correctly the system     
    1161                         z1(ig)=1./(za(ig,1)+zb(ig,1)+
    1162      &                         zb(ig,2)*(1.-zd(ig,2)))
    1163                         zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+
    1164      &                           zb(ig,2)*zc(ig,2) +
    1165      &                           (-zdqsdif_tot(ig)) *subtimestep)
    1166      &                           * z1(ig)
    1167                         zd(ig,1)=zb(ig,1)*z1(ig)
    1168                         zq1temp(ig)=zc(ig,1)+ zd(ig,1)
    1169      &                              *qsat_ssi(ig,islope)
    1170            
    1171 
    11721209
    11731210                    else
    … …  
    11791216     &                      zb(ig,2)*zc(ig,2) +
    11801217     &                      (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
    1181                             zq1temp(ig)=zc(ig,1)
     1218                         zq1temp(ig)=zc(ig,1)
    11821219                    endif ! Of subsurface <-> atmosphere exchange
    11831220                 endif ! sublimating more than available frost & surface - frost exchange
    … …  
    12141251     &                                 *(qsat(ig)-qsat_ssi(ig,islope))
    12151252
    1216 ! 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)               
     1253! 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)             
     1254
     1255                  if ((-zdqsdif_tot(ig)+zdqsdif_ssi_frost(ig,islope))
     1256     &            *subtimestep.ge.(zqsurf(ig))) then
     1257                 
     1258                     zdqsdif_ssi_frost(ig,islope) =
     1259     &                         zqsurf(ig)/subtimestep
     1260     &                         +   zdqsdif_tot(ig)
     1261                  endif
     1262                 
    12171263                    zdqsdif_tot(ig) = zdqsdif_tot(ig) -
    1218      &                        zdqsdif_ssi_frost(ig,islope)
     1264     &                        zdqsdif_ssi_frost(ig,islope)
     1265
    12191266                 endif ! watercaptag or frost at the surface
    12201267             endif ! interaction frost <-> subsurface ice
    12211268             
    12221269
     1270ccc Special case : schorgho study
     1271
     1272       if (h2o_ice_depth(ig,islope).gt.0) then
     1273                    call watersat(1,Tice(ig,islope),pplev(ig,1),
     1274     &                            qsat_ssi(ig,islope))
     1275                  qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
     1276     &                                       *qsat_ssi(ig,islope)
     1277           flux_schorgho(ig,islope) = d_coef(ig,islope)
     1278     &                               /h2o_ice_depth(ig,islope)
     1279     &                              *rho(ig)*dryness(ig)
     1280     &                              *(min(zq1temp(ig),qsat(ig))
     1281     &                               -qsat_ssi(ig,islope))
     1282     
     1283            if((zqsurf(ig).gt.tol_frost)) then
     1284            flux_schorgho_vfrost(ig,islope) = d_coef(ig,islope)
     1285     &                               /h2o_ice_depth(ig,islope)
     1286     &                              *rho(ig)*dryness(ig)
     1287     &                              *(qsat(ig)
     1288     &                               -qsat_ssi(ig,islope))
     1289             else
     1290            flux_schorgho_vfrost(ig,islope) = d_coef(ig,islope)
     1291     &                               /h2o_ice_depth(ig,islope)
     1292     &                              *rho(ig)*dryness(ig)
     1293     &                              *(zq1temp(ig)
     1294     &                               -qsat_ssi(ig,islope))
     1295                 
     1296             endif
     1297           
     1298        endif
     1299
     1300ccc
    12231301c             Starting upward calculations for water :
    12241302             zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
    … …  
    12461324     &                              zdqsdif_ssi_atm_tot(ig,islope)
    12471325     &                            + zdqsdif_ssi_atm(ig,islope)
     1326     &                            * subtimestep
    12481327              zdqsdif_ssi_frost_tot(ig,islope) =
    12491328     &                                zdqsdif_ssi_frost_tot(ig,islope)
    12501329     &                              + zdqsdif_ssi_frost(ig,islope)
     1330     &                              * subtimestep
     1331     
     1332              flux_schorgho_vfrost_tot(ig,islope) =
     1333     &           flux_schorgho_vfrost_tot(ig,islope) +
     1334     &           flux_schorgho_vfrost(ig,islope)* subtimestep
     1335
     1336              flux_schorgho_tot(ig,islope) =
     1337     &           flux_schorgho_tot(ig,islope) +
     1338     &           flux_schorgho(ig,islope)* subtimestep     
     1339     
     1340     
    12511341c             Monitoring instantaneous latent heat flux in W.m-2 :
    12521342              zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
    … …  
    12761366     &                       cos(pi*def_slope_mean(islope)/180.))
    12771367     &                       /ptimestep
    1278      
     1368   
     1369
     1370            zdqsoildif_tot(ig,:,:,islope) =  (zqsoil(ig,:,:,islope) -
     1371     &                      qsoil(ig,:,:,islope)/
     1372     &                       cos(pi*def_slope_mean(islope)/180.))
     1373     &                       /ptimestep
     1374            qsoil(ig,:,:,islope) = zqsoil(ig,:,:,islope)
    12791375            zdqsdif_ssi_tot(ig,islope) =
    12801376     &                        zdqsdif_ssi_atm_tot(ig,islope)
    … …  
    14131509     &                zdqsdif_ssi_tot(:,iflat))
    14141510
    1415 
     1511        call write_output('zdqsdif_regolith_frost_tot',
     1512     &                '','',zdqsdif_regolith_frost_tot(:,iflat))
     1513
     1514        call write_output('zdqsdif_regolith_atm_tot',
     1515     &                '','',zdqsdif_regolith_atm_tot(:,iflat))
     1516
     1517
     1518        call write_output('zdqsoildif_tot','','',
     1519     &           zdqsoildif_tot(:,:,igcm_h2o_ice_soil,iflat))
     1520     
     1521             call write_output('zdqsoildif_tot','','',
     1522     &           zdqsoildif_tot(:,:,igcm_h2o_ice_soil,iflat))
     1523   
     1524   
     1525         call write_output('flux_schorgho_vfrost_tot',
     1526     &     '','kg.m-2.s-1',
     1527     &                flux_schorgho_vfrost_tot(:,iflat))
     1528     
     1529         call write_output('flux_schorgho_tot',
     1530     &     '','kg.m-2.s-1',
     1531     &                flux_schorgho_tot(:,iflat))
    14161532C       Diagnostic output for HDO
    14171533!        if (hdo) then
    … …  
    14231539!     &                       ' ',h2oflux(:))
    14241540!        endif
     1541
     1542         call write_output("flux_soillayer","",
     1543     &             "",flux_soil_py_tot(:,:,iflat))
     1544
    14251545
    14261546c-----------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.