Ignore:
Timestamp:
Mar 12, 2024, 5:00:40 PM (10 months ago)
Author:
llange
Message:

Mars PCM
reverting the revert commit -r3260. I've commited the wrong trunk. Sorry.
The bug with the soil index is fixed.
LL

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

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

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

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

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

    r3261 r3262  
    4040     &                           lag_layer
    4141      use microphys_h, only: mteta
    42       use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D,
    43      & ads_massive_ice
     42      use comsoil_h, only: adsorption_soil, choice_ads
    4443      use nonoro_gwd_mix_mod, only: calljliu_gwimix
    4544
    … …  
    10431042         call getin_p("choice_ads",choice_ads)
    10441043         
    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)
    10511044         poreice_tifeedback = .false.
    10521045         call getin_p("poreice_tifeedback",poreice_tifeedback)

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

    r3261 r3262  
    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, ads_massive_ice
     24                                    ini_comsoil_h_slope_var, end_comsoil_h_slope_var
    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
    118 logical :: ice_fb = .false.
    119118
    120119!=======================================================================
    … …  
    670669if (.not. therestartfi) qsoil = 0.
    671670
    672 call getin("ice_fb_test1D",ice_fb)
    673671if (.not. therestartfi) then
    674672    ! Initialize depths
    … …  
    690688            if(adsorption_soil) then
    691689            ! add subsurface ice in qsoil if one runs with adsorption
    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
     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
    699692            endif
    700693        else ! searching for the ice/regolith boundary:
    701694            do isoil = 1,nsoil
    702695                if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
    703                     iref = isoil +1
     696                    iref = isoil + 1
    704697                    exit
    705698                endif
    … …  
    727720        inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
    728721    endif ! ice_depth > 0
    729 
    730 if(.not.(ice_fb)) then
    731 inertiedat(1,:) = inertiedat(1,1)
    732 end if
     722   
    733723    do isoil = 1,nsoil
    734724        inertiesoil(1,isoil,:) = inertiedat(1,isoil)

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

    r3261 r3262  
    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

    r3261 r3262  
    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, compute_Tice
     84      USE vdifc_mod, ONLY: vdifc
    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
    127126      IMPLICIT NONE
    128127c=======================================================================
    … …  
    543542! Variable for ice table
    544543      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]
    546544      REAL :: rhowater_surf_sat(ngrid,nslope)        ! Water density at the surface at saturation [kg/m^3]
    547545      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]
    549546      REAL,PARAMETER :: alpha_clap_h2o = 28.9074     ! Coeff for Clapeyron law [/]
    550547      REAL,PARAMETER :: beta_clap_h2o = -6143.7      ! Coeff for Clapeyron law [K]
    … …  
    554551      REAL :: ztmp1,ztmp2                            ! intermediate variables to compute the mean molar mass of the layer
    555552      REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores
    556       REAL :: Tice
    557553! Variable for the computation of the TKE with parameterization from ATKE working group
    558554      REAL :: viscom                              ! kinematic molecular viscosity for momentum
    … …  
    39433939     &         * mmol(igcm_h2o_vap)/(mugaz*r)
    39443940             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)
    39503941             DO isoil = 1,nsoilmx
    39513942             rhowater_soil(ig,isoil,islope) =
    … …  
    39543945     &         * mmol(igcm_h2o_vap)/(mugaz*r)
    39553946             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      
    39683947          ENDDO
    39693948        ENDDO
    … …  
    39753954     &     "rhowater_surface",'kg.m-3',
    39763955     &     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))
    39853956      DO islope = 1,nslope
    39863957        write(str2(1:2),'(i2.2)') islope

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

    r3261 r3262  
    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,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
     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
    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
    1312implicit none
    1413
    … …  
    184183logical, allocatable, save :: close_top(:, :)         ! closing diffusion at the top of a layer if ice rises over saturation
    185184logical, allocatable, save :: close_bottom(:, :)      ! closing diffusion at the bottom of a layer if ice risies over saturation
    186 logical, parameter :: print_closure_warnings = .false.
     185logical, parameter :: print_closure_warnings = .true.
    187186
    188187! Coefficients for the Diffusion calculations
    … …  
    270269real*8 :: diff(ngrid, nsoil)              ! difference between total water content and ice + vapor content
    271270                                          ! (only used for output, inconstistent: should just be adswater)
    272 real,intent(out) :: var_flux_soil(ngrid, nsoil)
     271real :: var_flux_soil(ngrid, nsoil)
    273272
    274273! Loop variables and counters
    … …  
    569568                  ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 -  only exact for normal diffusion)
    570569                  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))
    571                   if(ads_const_D) D_mid(ig, ik) = d_coef(ig,1)
     570
    572571                  ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991)
    573572                  !                             D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik)  &
    … …  
    594593                 
    595594                  saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik))  ! new diffusion interaction
    596                  
     595                  
    597596                  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                  
    601597                 
    602598!                  D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) &  ! old implementation with direct interpolation
    … …  
    15711567!       call write_output("dztot1","change in dztot", "unknown",real(dztot1))
    15721568
     1569        call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil)                               
    15731570     
    15741571!       call write_output("A_coef","A coefficient of dztot1", "unknown",real(B))                               
    … …  
    15981595       call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater))                   
    15991596     
    1600         call write_output("subsurfice_mass_py","subsurface ice", "kg/m^3",real(ice))                   
     1597!       call write_output("subsurfice","subsurface ice", "kg/m^3",real(ice))                   
    16011598
    16021599        call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego)                   
    … …  
    16061603!       call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out))                         
    16071604
    1608        call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_mid)                         
     1605!       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter)                         
    16091606           
    16101607!endif

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

    r3261 r3262  
    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,
    38      &                     igcm_h2o_ice_soil
     37      use comsoil_h, only: layer, mlayer,adsorption_soil
    3938      use vdif_cd_mod, only: vdif_cd
    4039      use lmdz_call_atke, only: call_atke
    … …  
    235234      LOGICAL :: writeoutput          ! boolean to say to soilexchange.F if we are at the last iteration and thus if he  can write in the diagsoil
    236235      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)
    240236!! Water buyoncy
    241237      LOGICAL :: virtual
    … …  
    250246      REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep
    251247      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)     
    262248     
    263249      REAL :: tol_frost = 1e-4 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect
    … …  
    320306      dwatercap_dif(1:ngrid,1:nslope)=0
    321307      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
    324308      zq1temp_regolith(1:ngrid)=0
    325309      zdqsdif_tot(1:ngrid)=0
    … …  
    331315      zdqsdif_ssi_atm_tot(:,:) = 0.
    332316      zdqsdif_ssi_frost_tot(:,:) = 0.
    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.
    341317     
    342      
     318           
    343319c    ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp
    344320c       avec rho=p/RT=p/ (R Theta) (p/ps)**kappa
    … …  
    10281004              watercap_tmp(ig) = watercap(ig,islope)/
    10291005     &                       cos(pi*def_slope_mean(islope)/180.)
    1030 
    1031               zqsoil(ig,:,:,islope) = qsoil(ig,:,:,islope)/
    1032      &             cos(pi*def_slope_mean(islope)/180.)
    1033 
    10341006           ENDDO ! ig=1,ngrid
    10351007           
    … …  
    10431015           saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap)   
    10441016           DO ig=1,ngrid
    1045 !           nsubtimestep(ig)=1 !for debug
     1017c           nsubtimestep(ig)=1 !for debug
    10461018           subtimestep = ptimestep/nsubtimestep(ig)
    10471019                          call write_output('subtimestep',
    … …  
    11171089     &                     za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:),
    11181090     &                     zdqsdif_surf(ig), zqsurf(ig),
    1119      &                     zqsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
     1091     &                     qsoil(ig,:,:,islope), pplev(ig,1), rho(ig),
    11201092     &                     writeoutput,zdqsdif_regolith(ig,islope),
    11211093     &                     zq1temp_regolith(ig),
    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
     1094     &                     pore_icefraction(ig,:,islope))
     1095
    11271096
    11281097                 if(.not.watercaptag(ig)) then
    … …  
    11311100                         zdqsdif_tot(ig)=
    11321101     &                        -zqsurf(ig)/subtimestep
    1133                          zdqsdif_regolith_atm_tot(ig,islope) =
    1134      &                         zdqsdif_regolith_atm_tot(ig,islope)
    1135      &                   + zdqsdif_regolith(ig,islope)*subtimestep
    11361102                     else
    11371103                         zdqsdif_tot(ig) = zdqsdif_surf(ig) +
    11381104     &                           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
    11421105                     endif  ! of "if exchange = true"
    11431106                 endif  ! of "if not.watercaptag"
    … …  
    11831146                         qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope)
    11841147     &                                       *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            
    11971148                  ! Flux from the subsurface     
    11981149                        if(old_wsublimation_scheme) then
    … …  
    12071158     &                     *(zq1temp(ig)-qsat_ssi(ig,islope))                         
    12081159                        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
    12091172
    12101173                    else
    … …  
    12161179     &                      zb(ig,2)*zc(ig,2) +
    12171180     &                      (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)
    1218                          zq1temp(ig)=zc(ig,1)
     1181                            zq1temp(ig)=zc(ig,1)
    12191182                    endif ! Of subsurface <-> atmosphere exchange
    12201183                 endif ! sublimating more than available frost & surface - frost exchange
    … …  
    12511214     &                                 *(qsat(ig)-qsat_ssi(ig,islope))
    12521215
    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                  
     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)               
    12631217                    zdqsdif_tot(ig) = zdqsdif_tot(ig) -
    1264      &                        zdqsdif_ssi_frost(ig,islope)
    1265 
     1218     &                        zdqsdif_ssi_frost(ig,islope)
    12661219                 endif ! watercaptag or frost at the surface
    12671220             endif ! interaction frost <-> subsurface ice
    12681221             
    12691222
    1270 ccc 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 
    1300 ccc
    13011223c             Starting upward calculations for water :
    13021224             zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig)
    … …  
    13241246     &                              zdqsdif_ssi_atm_tot(ig,islope)
    13251247     &                            + zdqsdif_ssi_atm(ig,islope)
    1326      &                            * subtimestep
    13271248              zdqsdif_ssi_frost_tot(ig,islope) =
    13281249     &                                zdqsdif_ssi_frost_tot(ig,islope)
    13291250     &                              + 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      
    13411251c             Monitoring instantaneous latent heat flux in W.m-2 :
    13421252              zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+
    … …  
    13661276     &                       cos(pi*def_slope_mean(islope)/180.))
    13671277     &                       /ptimestep
    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)
     1278     
    13751279            zdqsdif_ssi_tot(ig,islope) =
    13761280     &                        zdqsdif_ssi_atm_tot(ig,islope)
    … …  
    15091413     &                zdqsdif_ssi_tot(:,iflat))
    15101414
    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))
     1415
    15321416C       Diagnostic output for HDO
    15331417!        if (hdo) then
    … …  
    15391423!     &                       ' ',h2oflux(:))
    15401424!        endif
    1541 
    1542          call write_output("flux_soillayer","",
    1543      &             "",flux_soil_py_tot(:,:,iflat))
    1544 
    15451425
    15461426c-----------------------------------------------------------------------

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

    r3250 r3262  
    107107              DO ik=1,nsoil-1
    108108                 IF ((icedepth.ge.layer(ik)).and.icedepth.lt.layer(ik+1)) THEN
    109                     iref=ik
     109                    iref=ik+1
    110110                 EXIT
    111111                 ENDIF
    … …  
    116116              ENDDO
    117117!             Transition (based on the equations of thermal conduction):
    118               newtherm_i(ig,iref,islope)=sqrt( (layer(iref+1)-layer(iref))/(((icedepth-layer(iref))/inert_h2o_ice**2) + &
    119               ((layer(iref+1)-icedepth)/inertiedat(ig,ik)**2) ) )
     118              newtherm_i(ig,iref,islope)=sqrt( (layer(iref)-layer(iref-1))/(((icedepth-layer(iref-1))/inert_h2o_ice**2) + &
     119              ((layer(iref)-icedepth)/inertiedat(ig,ik)**2) ) )
    120120!              Underlying regolith:
    121121              DO ik=iref+1,nsoil
Note: See TracChangeset for help on using the changeset viewer.