Index: trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F (revision 3261) @@ -5,5 +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 Index: trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90 (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90 (revision 3261) @@ -31,4 +31,6 @@ ! 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,4 +41,5 @@ REAL, parameter :: porosity_reg = 0.45 !$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads) +!$OMP THREADPRIVATE(ads_const_D,ads_massive_ice) !======================================================================= Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 3261) @@ -40,5 +40,6 @@ & 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,4 +1043,10 @@ 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) Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90 (revision 3261) @@ -22,5 +22,5 @@ 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,4 +116,5 @@ ! LL: Subsurface water ice real :: rho_H2O_ice ! Density (kg/m^3) of water ice +logical :: ice_fb = .false. !======================================================================= @@ -669,4 +670,5 @@ if (.not. therestartfi) qsoil = 0. +call getin("ice_fb_test1D",ice_fb) if (.not. therestartfi) then ! Initialize depths @@ -688,11 +690,16 @@ 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,5 +708,5 @@ 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,6 +715,11 @@ ! 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,5 +727,8 @@ 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) Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90 (revision 3261) @@ -257,5 +257,5 @@ ! 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 Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (revision 3261) @@ -82,5 +82,5 @@ 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,4 +124,5 @@ 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,6 +543,8 @@ ! 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,4 +554,5 @@ 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,4 +3943,9 @@ & * 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,4 +3954,16 @@ & * 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,4 +3975,12 @@ & "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 Index: trunk/LMDZ.MARS/libf/phymars/soilwater.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/soilwater.F90 (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/soilwater.F90 (revision 3261) @@ -1,8 +1,8 @@ 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,4 +10,5 @@ use geometry_mod, only: cell_area, latitude_deg use write_output_mod, only: write_output + use paleoclimate_mod, only: d_coef implicit none @@ -183,5 +184,5 @@ 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,5 +270,5 @@ 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,5 +569,5 @@ ! 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,6 +594,9 @@ 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,5 +1571,4 @@ ! 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,5 +1598,5 @@ 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,5 +1606,5 @@ ! 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 Index: trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (revision 3253) +++ trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F (revision 3261) @@ -35,5 +35,6 @@ 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,4 +235,7 @@ 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,4 +250,14 @@ 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,4 +320,6 @@ 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,6 +331,14 @@ 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,4 +1028,8 @@ 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,5 +1043,5 @@ 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,9 +1117,12 @@ & 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,7 +1131,13 @@ 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,4 +1183,16 @@ 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,16 +1207,4 @@ & *(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,5 +1216,5 @@ & 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,11 +1251,52 @@ & *(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,7 +1324,19 @@ & 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,5 +1366,11 @@ & 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,5 +1509,25 @@ & 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,4 +1539,8 @@ ! & ' ',h2oflux(:)) ! endif + + call write_output("flux_soillayer","", + & "",flux_soil_py_tot(:,:,iflat)) + c-----------------------------------------------------------------------