- Timestamp:
- Mar 12, 2024, 5:00:40 PM (8 months ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 9 edited
-
co2condens_mod.F (modified) (1 diff)
-
comsoil_h.F90 (modified) (2 diffs)
-
conf_phys.F (modified) (2 diffs)
-
dyn1d/init_testphys1d_mod.F90 (modified) (5 diffs)
-
dyn1d/testphys1d.F90 (modified) (1 diff)
-
physiq_mod.F (modified) (7 diffs)
-
soilwater.F90 (modified) (9 diffs)
-
vdifc_mod.F (modified) (17 diffs)
-
waterice_tifeedback_mod.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
r3261 r3262 5 5 logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2 6 6 !$OMP THREADPRIVATE(scavco2cond) 7 real, save :: CO2cond_ps = 0. ! Coefficient to control the surface pressure change7 real, save :: CO2cond_ps = 1. ! Coefficient to control the surface pressure change 8 8 9 9 CONTAINS -
trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
r3261 r3262 31 31 ! Subsurface tracers: 32 32 logical, save :: adsorption_soil ! boolean to call adosrption (or not) 33 logical, save :: ads_const_D ! boolean to have constant diffusion coefficient34 logical, save :: ads_massive_ice ! boolean to have massive subsurface ice35 33 real, save :: choice_ads ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90) 36 34 real, save, allocatable, dimension(:,:,:,:) :: qsoil ! subsurface tracers (kg/m^3 of regol) … … 41 39 REAL, parameter :: porosity_reg = 0.45 42 40 !$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads) 43 !$OMP THREADPRIVATE(ads_const_D,ads_massive_ice)44 41 45 42 !======================================================================= -
trunk/LMDZ.MARS/libf/phymars/conf_phys.F
r3261 r3262 40 40 & lag_layer 41 41 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 44 43 use nonoro_gwd_mix_mod, only: calljliu_gwimix 45 44 … … 1043 1042 call getin_p("choice_ads",choice_ads) 1044 1043 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)1051 1044 poreice_tifeedback = .false. 1052 1045 call getin_p("poreice_tifeedback",poreice_tifeedback) -
trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
r3261 r3262 22 22 use comsoil_h, only: volcapa, nsoilmx, inertiesoil, inertiedat, layer, mlayer, flux_geo, tsoil, & 23 23 adsorption_soil, qsoil, igcm_h2o_ice_soil, porosity_reg, & 24 ini_comsoil_h_slope_var, end_comsoil_h_slope_var , ads_massive_ice24 ini_comsoil_h_slope_var, end_comsoil_h_slope_var 25 25 use comvert_mod, only: ap, bp, aps, bps, pa, preff, presnivs, pseudoalt, scaleheight 26 26 use dimradmars_mod, only: tauvis, totcloudfrac, albedo, ini_dimradmars_mod_slope_var, end_dimradmars_mod_slope_var … … 116 116 ! LL: Subsurface water ice 117 117 real :: rho_H2O_ice ! Density (kg/m^3) of water ice 118 logical :: ice_fb = .false.119 118 120 119 !======================================================================= … … 670 669 if (.not. therestartfi) qsoil = 0. 671 670 672 call getin("ice_fb_test1D",ice_fb)673 671 if (.not. therestartfi) then 674 672 ! Initialize depths … … 690 688 if(adsorption_soil) then 691 689 ! 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 699 692 endif 700 693 else ! searching for the ice/regolith boundary: 701 694 do isoil = 1,nsoil 702 695 if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then 703 iref = isoil + 1696 iref = isoil + 1 704 697 exit 705 698 endif … … 727 720 inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia 728 721 endif ! ice_depth > 0 729 730 if(.not.(ice_fb)) then 731 inertiedat(1,:) = inertiedat(1,1) 732 end if 722 733 723 do isoil = 1,nsoil 734 724 inertiesoil(1,isoil,:) = inertiedat(1,isoil) -
trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
r3261 r3262 257 257 ! Compute pressure for next time step 258 258 ! ----------------------------------- 259 psurf = psurf !+ dttestphys*dpsurf(1) ! Surface pressure change259 psurf = psurf + dttestphys*dpsurf(1) ! Surface pressure change 260 260 plev = ap + psurf*bp 261 261 play = aps + psurf*bps -
trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
r3261 r3262 82 82 USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 83 83 USE calldrag_noro_mod, ONLY: calldrag_noro 84 USE vdifc_mod, ONLY: vdifc , compute_Tice84 USE vdifc_mod, ONLY: vdifc 85 85 use param_v4_h, only: nreact,n_avog, 86 86 & fill_data_thermos, allocate_param_thermos … … 124 124 use lmdz_atke_turbulence_ini, only : atke_ini 125 125 use waterice_tifeedback_mod, only : waterice_tifeedback 126 use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer127 126 IMPLICIT NONE 128 127 c======================================================================= … … 543 542 ! Variable for ice table 544 543 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]546 544 REAL :: rhowater_surf_sat(ngrid,nslope) ! Water density at the surface at saturation [kg/m^3] 547 545 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]549 546 REAL,PARAMETER :: alpha_clap_h2o = 28.9074 ! Coeff for Clapeyron law [/] 550 547 REAL,PARAMETER :: beta_clap_h2o = -6143.7 ! Coeff for Clapeyron law [K] … … 554 551 REAL :: ztmp1,ztmp2 ! intermediate variables to compute the mean molar mass of the layer 555 552 REAL :: pore_icefraction(ngrid,nsoilmx,nslope) ! ice filling fraction in the pores 556 REAL :: Tice557 553 ! Variable for the computation of the TKE with parameterization from ATKE working group 558 554 REAL :: viscom ! kinematic molecular viscosity for momentum … … 3943 3939 & * mmol(igcm_h2o_vap)/(mugaz*r) 3944 3940 endif 3945 ! schorghofer way3946 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)3950 3941 DO isoil = 1,nsoilmx 3951 3942 rhowater_soil(ig,isoil,islope) = … … 3954 3945 & * mmol(igcm_h2o_vap)/(mugaz*r) 3955 3946 ENDDO 3956 if(h2o_ice_depth(ig,islope).gt.0) then3957 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 & / Tice3963 & * mmol(igcm_h2o_vap)/(mugaz*r)3964 else3965 rhowater_ice(ig,islope) = 0.3966 endif3967 3968 3947 ENDDO 3969 3948 ENDDO … … 3975 3954 & "rhowater_surface",'kg.m-3', 3976 3955 & 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))3985 3956 DO islope = 1,nslope 3986 3957 write(str2(1:2),'(i2.2)') islope -
trunk/LMDZ.MARS/libf/phymars/soilwater.F90
r3261 r3262 1 1 subroutine soilwater(ngrid, nlayer, nq, nsoil, nqsoil, ptsrf, ptsoil, ptimestep, & 2 2 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_D3 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 7 7 use comcstfi_h 8 8 use tracer_mod … … 10 10 use geometry_mod, only: cell_area, latitude_deg 11 11 use write_output_mod, only: write_output 12 use paleoclimate_mod, only: d_coef13 12 implicit none 14 13 … … 184 183 logical, allocatable, save :: close_top(:, :) ! closing diffusion at the top of a layer if ice rises over saturation 185 184 logical, allocatable, save :: close_bottom(:, :) ! closing diffusion at the bottom of a layer if ice risies over saturation 186 logical, parameter :: print_closure_warnings = . false.185 logical, parameter :: print_closure_warnings = .true. 187 186 188 187 ! Coefficients for the Diffusion calculations … … 270 269 real*8 :: diff(ngrid, nsoil) ! difference between total water content and ice + vapor content 271 270 ! (only used for output, inconstistent: should just be adswater) 272 real ,intent(out):: var_flux_soil(ngrid, nsoil)271 real :: var_flux_soil(ngrid, nsoil) 273 272 274 273 ! Loop variables and counters … … 569 568 ! calculate the midlayer diffusion coeff. (with dependence on saturation_water_ice from Meslin et al., 2010 - only exact for normal diffusion) 570 569 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 572 571 ! Midlayer diffusion coeff. (correlation by Rogers and Nielson, 1991) 573 572 ! D_mid(ig, ik) = D0(ig, ik) * (1. - saturation_water_ice(ig, ik)) * exp(-6. * saturation_water_ice(ig, ik) & … … 594 593 595 594 saturation_water_ice_inter(ig, ik) = min(saturation_water_ice(ig, ik + 1), saturation_water_ice(ig, ik)) ! new diffusion interaction 596 595 597 596 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 601 597 602 598 ! D_inter(ig, ik) = ( (layer(ik) - mlayer(ik - 1)) * D_mid(ig, ik + 1) & ! old implementation with direct interpolation … … 1571 1567 ! call write_output("dztot1","change in dztot", "unknown",real(dztot1)) 1572 1568 1569 call write_output("flux_soillayer","flux of water between the soil layers", "kg m-2 s-1",var_flux_soil) 1573 1570 1574 1571 ! call write_output("A_coef","A coefficient of dztot1", "unknown",real(B)) … … 1598 1595 call write_output("adswater","subsurface adsorbed water", "kg/m^3",real(adswater)) 1599 1596 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)) 1601 1598 1602 1599 call write_output("flux_rego",'flux of water from the regolith','kg/m^2',zdqsdifrego) … … 1606 1603 ! call write_output("close",'close top, bottom, or none (1, -1, 0)','no unit',real(close_out)) 1607 1604 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) 1609 1606 1610 1607 !endif -
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r3261 r3262 35 35 use microphys_h, only: To 36 36 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 39 38 use vdif_cd_mod, only: vdif_cd 40 39 use lmdz_call_atke, only: call_atke … … 235 234 LOGICAL :: writeoutput ! boolean to say to soilexchange.F if we are at the last iteration and thus if he can write in the diagsoil 236 235 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)240 236 !! Water buyoncy 241 237 LOGICAL :: virtual … … 250 246 REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep 251 247 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)262 248 263 249 REAL :: tol_frost = 1e-4 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect … … 320 306 dwatercap_dif(1:ngrid,1:nslope)=0 321 307 zdqsdif_regolith(1:ngrid,1:nslope)=0 322 zdqsdif_regolith_frost_tot(1:ngrid,1:nslope)=0323 zdqsdif_regolith_atm_tot(1:ngrid,1:nslope)=0324 308 zq1temp_regolith(1:ngrid)=0 325 309 zdqsdif_tot(1:ngrid)=0 … … 331 315 zdqsdif_ssi_atm_tot(:,:) = 0. 332 316 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.341 317 342 318 343 319 c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp 344 320 c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa … … 1028 1004 watercap_tmp(ig) = watercap(ig,islope)/ 1029 1005 & cos(pi*def_slope_mean(islope)/180.) 1030 1031 zqsoil(ig,:,:,islope) = qsoil(ig,:,:,islope)/1032 & cos(pi*def_slope_mean(islope)/180.)1033 1034 1006 ENDDO ! ig=1,ngrid 1035 1007 … … 1043 1015 saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap) 1044 1016 DO ig=1,ngrid 1045 !nsubtimestep(ig)=1 !for debug1017 c nsubtimestep(ig)=1 !for debug 1046 1018 subtimestep = ptimestep/nsubtimestep(ig) 1047 1019 call write_output('subtimestep', … … 1117 1089 & za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:), 1118 1090 & zdqsdif_surf(ig), zqsurf(ig), 1119 & zqsoil(ig,:,:,islope), pplev(ig,1), rho(ig),1091 & qsoil(ig,:,:,islope), pplev(ig,1), rho(ig), 1120 1092 & writeoutput,zdqsdif_regolith(ig,islope), 1121 1093 & 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 1127 1096 1128 1097 if(.not.watercaptag(ig)) then … … 1131 1100 zdqsdif_tot(ig)= 1132 1101 & -zqsurf(ig)/subtimestep 1133 zdqsdif_regolith_atm_tot(ig,islope) =1134 & zdqsdif_regolith_atm_tot(ig,islope)1135 & + zdqsdif_regolith(ig,islope)*subtimestep1136 1102 else 1137 1103 zdqsdif_tot(ig) = zdqsdif_surf(ig) + 1138 1104 & 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)*subtimestep1142 1105 endif ! of "if exchange = true" 1143 1106 endif ! of "if not.watercaptag" … … 1183 1146 qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope) 1184 1147 & *qsat_ssi(ig,islope) 1185 1186 ! And now we solve correctly the system1187 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 1197 1148 ! Flux from the subsurface 1198 1149 if(old_wsublimation_scheme) then … … 1207 1158 & *(zq1temp(ig)-qsat_ssi(ig,islope)) 1208 1159 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 1209 1172 1210 1173 else … … 1216 1179 & zb(ig,2)*zc(ig,2) + 1217 1180 & (-zdqsdif_tot(ig)) *subtimestep) *z1(ig) 1218 zq1temp(ig)=zc(ig,1)1181 zq1temp(ig)=zc(ig,1) 1219 1182 endif ! Of subsurface <-> atmosphere exchange 1220 1183 endif ! sublimating more than available frost & surface - frost exchange … … 1251 1214 & *(qsat(ig)-qsat_ssi(ig,islope)) 1252 1215 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) 1263 1217 zdqsdif_tot(ig) = zdqsdif_tot(ig) - 1264 & zdqsdif_ssi_frost(ig,islope) 1265 1218 & zdqsdif_ssi_frost(ig,islope) 1266 1219 endif ! watercaptag or frost at the surface 1267 1220 endif ! interaction frost <-> subsurface ice 1268 1221 1269 1222 1270 ccc Special case : schorgho study1271 1272 if (h2o_ice_depth(ig,islope).gt.0) then1273 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)) then1284 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 else1290 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 endif1297 1298 endif1299 1300 ccc1301 1223 c Starting upward calculations for water : 1302 1224 zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig) … … 1324 1246 & zdqsdif_ssi_atm_tot(ig,islope) 1325 1247 & + zdqsdif_ssi_atm(ig,islope) 1326 & * subtimestep1327 1248 zdqsdif_ssi_frost_tot(ig,islope) = 1328 1249 & zdqsdif_ssi_frost_tot(ig,islope) 1329 1250 & + zdqsdif_ssi_frost(ig,islope) 1330 & * subtimestep1331 1332 flux_schorgho_vfrost_tot(ig,islope) =1333 & flux_schorgho_vfrost_tot(ig,islope) +1334 & flux_schorgho_vfrost(ig,islope)* subtimestep1335 1336 flux_schorgho_tot(ig,islope) =1337 & flux_schorgho_tot(ig,islope) +1338 & flux_schorgho(ig,islope)* subtimestep1339 1340 1341 1251 c Monitoring instantaneous latent heat flux in W.m-2 : 1342 1252 zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+ … … 1366 1276 & cos(pi*def_slope_mean(islope)/180.)) 1367 1277 & /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 1375 1279 zdqsdif_ssi_tot(ig,islope) = 1376 1280 & zdqsdif_ssi_atm_tot(ig,islope) … … 1509 1413 & zdqsdif_ssi_tot(:,iflat)) 1510 1414 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 1532 1416 C Diagnostic output for HDO 1533 1417 ! if (hdo) then … … 1539 1423 ! & ' ',h2oflux(:)) 1540 1424 ! endif 1541 1542 call write_output("flux_soillayer","",1543 & "",flux_soil_py_tot(:,:,iflat))1544 1545 1425 1546 1426 c----------------------------------------------------------------------- -
trunk/LMDZ.MARS/libf/phymars/waterice_tifeedback_mod.F90
r3250 r3262 107 107 DO ik=1,nsoil-1 108 108 IF ((icedepth.ge.layer(ik)).and.icedepth.lt.layer(ik+1)) THEN 109 iref=ik 109 iref=ik+1 110 110 EXIT 111 111 ENDIF … … 116 116 ENDDO 117 117 ! 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) ) ) 120 120 ! Underlying regolith: 121 121 DO ik=iref+1,nsoil
Note: See TracChangeset
for help on using the changeset viewer.
