Changeset 3253 for trunk/LMDZ.MARS
- Timestamp:
- Mar 5, 2024, 11:36:32 AM (9 months ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 4 edited
-
changelog.txt (modified) (1 diff)
-
deftank/field_def_physics_mars.xml (modified) (2 diffs)
-
libf/phymars/dyn1d/init_testphys1d_mod.F90 (modified) (4 diffs)
-
libf/phymars/vdifc_mod.F (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/changelog.txt
r3251 r3253 4541 4541 == 04/03/2024 == LL 4542 4542 Introduce the possibility to prescribe (in 1D) the soil water ice content when running with adsorption/desoprtion 4543 4544 == 05/03/2024 == LL 4545 Following -r 3098; cleaning of vdfic for the management of subsurface water ice. 4546 Fixing some errors (wrong interpolation to compute the water ice temperature, wrong boundary conditions to compute qvap(1)) -
trunk/LMDZ.MARS/deftank/field_def_physics_mars.xml
r3136 r3253 533 533 long_name="co2 in the first layer" 534 534 unit="kg/kg" /> 535 536 535 <!-- Subsurface ice interactions --> 536 <field id="zdqsdif_ssi_frost_tot" 537 long_name="Flux between frost and subsurface ice" 538 unit="kg.m-2.s-1" /> 539 <field id="zdqsdif_ssi_atm_tot" 540 long_name="Flux between atmosphere and subsurface ice" 541 unit="kg.m-2.s-1" /> 542 <field id="zdqsdif_ssi_tot" 543 long_name="Total flux echanged with subsurface ice" 544 unit="kg.m-2.s-1" /> 537 545 538 546 <!-- CO2 cycle --> … … 1255 1263 <field id="waterdensity_soil_slope07" 1256 1264 long_name="waterdensity_soil of slope 07" 1257 unit="kg.m-3" /> 1258 <field id="zdqsdif_ssi_frost" 1259 long_name="Flux between frost and subsurface" 1260 unit="kg.m-2.s-1" /> 1265 unit="kg.m-3" /> 1261 1266 <field id="subtimestep" 1262 1267 long_name="vdifc substimestep length" -
trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
r3251 r3253 42 42 use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd 43 43 use conf_phys_mod, only: conf_phys 44 use paleoclimate_mod, only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h 44 use paleoclimate_mod, only: paleoclimate, ini_paleoclimate_h, end_paleoclimate_h, h2o_ice_depth 45 45 use comslope_mod, only: nslope, subslope_dist, ini_comslope_h, end_comslope_h 46 46 use co2condens_mod, only: CO2cond_ps … … 541 541 call getin("subsurface_ice_depth",ice_depth) 542 542 write(*,*) " subsurface_ice_depth = ",ice_depth 543 543 h2o_ice_depth(:,:) = ice_depth 544 544 z0(1) = z0_default ! default value for roughness 545 545 write(*,*) 'Surface roughness length z0 (m)?' … … 668 668 669 669 if (.not. therestartfi) qsoil = 0. 670 671 670 672 671 if (.not. therestartfi) then … … 716 715 inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia 717 716 endif ! ice_depth > 0 718 717 719 718 do isoil = 1,nsoil 720 719 inertiesoil(1,isoil,:) = inertiedat(1,isoil) -
trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F
r3248 r3253 128 128 REAL zcst1 129 129 REAL zu2(ngrid) 130 REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located. 131 REAL qeq(ngrid,nslope) ! saturation water vapor in the subsurface 132 REAL dist_up(ngrid,nslope) !distance from ice to layer above 133 REAL dist_down(ngrid,nslope) !distance from ice to layer down 134 REAL dist_sum(ngrid,nslope) ! sum of distance 135 REAL zdqsdif_ssi(ngrid,nslope) !SSI flux 136 REAL zdqsdif_ssi_frost(ngrid,nslope) !SSI-frost interaction 130 137 131 138 132 EXTERNAL SSUM,SCOPY … … 172 166 REAL rho(ngrid) ! near surface air density 173 167 REAL qsat(ngrid) 174 REAL qsat2(ngrid,nslope) 175 REAL resist(ngrid,nslope) !subsurface ice flux reduction coef 168 176 169 177 170 REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO … … 244 237 LOGICAL :: virtual 245 238 246 247 239 !! Subsurface ice interactions 240 REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located (K) 241 REAL qsat_ssi(ngrid,nslope) ! qsat(Tice) (kg/kg) 242 REAL resist(ngrid,nslope) ! subsurface ice flux reduction coef (1) 243 REAL zdqsdif_ssi_atm(ngrid,nslope) ! SSI - atmosphere flux (kg/m^2/s^-1) at a given subtimestep 244 REAL zdqsdif_ssi_frost(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) at a given subtimestep 245 REAL zdqsdif_ssi_atm_tot(ngrid,nslope) ! SSI - atmosphere flux (kg/m^2/s^-1) summed over all subtimestep 246 REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep 247 REAL zdqsdif_ssi_tot(ngrid,nslope) ! Total flux of the interactions with SSI kg/m^2/s^-1) 248 249 REAL :: tol_frost = 1e-4 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect 248 250 c ** un petit test de coherence 249 251 c -------------------------- … … 306 308 zq1temp_regolith(1:ngrid)=0 307 309 zdqsdif_tot(1:ngrid)=0 308 h2o_ice_depth(1:ngrid,1:nslope)=1309 310 virtual = .false. 310 311 pore_icefraction(:,:,:) = 0. 312 zdqsdif_ssi_atm(:,:) = 0. 313 zdqsdif_ssi_frost(:,:) = 0. 314 zdqsdif_ssi_tot(:,:) = 0. 315 zdqsdif_ssi_atm_tot(:,:) = 0. 316 zdqsdif_ssi_frost_tot(:,:) = 0. 317 318 311 319 c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp 312 320 c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa … … 983 991 c --------- h2o_vap -------------------------------- 984 992 985 986 c Traitement de la vapeur d'eau h2o_vap 987 c Utilisation d'un sous pas de temps afin 988 c de decrire le flux de chaleur latente 993 c Treatement of the water frost 994 c We use a subtimestep to take into account the release of latent heat 989 995 990 996 do iq=1,nq … … 993 999 DO islope = 1,nslope 994 1000 DO ig=1,ngrid 995 996 997 1001 998 1002 zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/ … … 1001 1005 & cos(pi*def_slope_mean(islope)/180.) 1002 1006 ENDDO ! ig=1,ngrid 1003 1004 c make_tsub : sous pas de temps adaptatif 1005 c la subroutine est a la fin du fichier 1007 1008 c Computation of the subtimestep 1006 1009 call make_tsub(ngrid,pdtsrf(:,islope),zqsurf, 1007 1010 & ptimestep,dtmax,watercaptag, 1008 & nsubtimestep) 1009 c Calculation for turbulent exchange with the surface (for ice)1011 & nsubtimestep) 1012 c Calculation for turbulent exchange (rho Cd,h U (qatm - qsat(Tsurf)) with the surface (for ice) 1010 1013 c initialization of ztsrf, which is surface temperature in 1011 1014 c the subtimestep. 1012 1015 saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap) 1013 1016 DO ig=1,ngrid 1014 !nsubtimestep(ig)=1 !for debug1015 1017 c nsubtimestep(ig)=1 !for debug 1018 subtimestep = ptimestep/nsubtimestep(ig) 1016 1019 call write_output('subtimestep', 1017 1020 & 'vdifc substimestep length','s',subtimestep) 1018 1021 ztsrf(ig)=ptsrf(ig,islope) ! +pdtsrf(ig)*subtimestep 1019 1022 zq_tmp_vap(ig,:,:) =zq(ig,:,:) 1020 c Debut du sous pas de temps1023 c Beginning of the subtimestep 1021 1024 DO tsub=1,nsubtimestep(ig) 1022 if(tsub.eq.nsubtimestep(ig)) writeoutput = .true. 1023 c C'est parti ! 1025 if(tsub.eq.nsubtimestep(ig)) writeoutput = .true. 1024 1026 zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) 1025 1027 & /float(nsubtimestep(ig)) … … 1061 1063 1062 1064 zdqsdif_tot(ig) = zdqsdif_surf(ig) 1063 !!! Subsurface exchange 1064 ! Check for subsurface exchanges 1065 ! -------------------------------------------------------------------------------------------------------------------------------- 1066 ! We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost 1067 ! when computing the complete adsorption/desorption model 1068 ! -------------------------------------------------------------------------------------------------------------------------------- 1065 1069 if(.not.watercaptag(ig)) then 1066 1070 if (((-(zdqsdif_surf(ig))* … … 1074 1078 exchange = .false. 1075 1079 endif 1076 zdqsdif_tot(ig) = zdqsdif_surf(ig) 1077 1080 1081 ! -------------------------------------------------------------------------------------------------------------------------------- 1082 ! If one consider adsorption, all the fluxes to/from the surface/subsurface/atmosphere are computed here 1083 ! -------------------------------------------------------------------------------------------------------------------------------- 1078 1084 1079 1085 if (adsorption_soil) then … … 1087 1093 & zq1temp_regolith(ig), 1088 1094 & pore_icefraction(ig,:,islope)) 1089 1090 1095 1091 1096 … … 1101 1106 endif ! of "if not.watercaptag" 1102 1107 endif ! adsorption 1103 1104 if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then 1105 if ((-zdqsdif_tot(ig)*subtimestep) 1106 & .gt.(zqsurf(ig))) then 1107 1108 ! -------------------------------------------------------------------------------------------------------------------------------------- 1109 ! Here we do the same, but without computing the complete adsorpption/desorption. Note that it work only if one does not use adsorption 1110 ! If no subsurface ice, then the models computes the surface flux/water vapor flux as usual 1111 ! -------------------------------------------------------------------------------------------------------------------------------------- 1112 1113 ! ------------------------------------------------------------------------------------------------------------------------------------------ 1114 ! First, We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost 1115 ! ------------------------------------------------------------------------------------------------------------------------------------------ 1116 1117 if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then 1118 if ((-zdqsdif_tot(ig)*subtimestep) 1119 & .ge.(zqsurf(ig))) then 1120 1121 zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep 1122 ! zqsurf(ig) = 0 1123 if((h2o_ice_depth(ig,islope).gt.0) 1124 & .and.lag_layer) then 1125 ! Atm <-> subsurface exchange, we need to update the exchange coefficient zb by a factor 1/1+R; R = zice*Cd,h*/D and add the flux from the subsurface 1126 1127 if(old_wsublimation_scheme) then 1128 resist(ig,islope)=h2o_ice_depth(ig,islope) 1129 & *zcdv(ig,islope)/d_coef(ig,islope) 1130 else 1131 resist(ig,islope)=h2o_ice_depth(ig,islope) 1132 & *zcdh(ig,islope)/d_coef(ig,islope) 1133 endif 1134 1135 zb(ig,1)=zb(ig,1)*1/(1+resist(ig,islope)) ! change zb to account subsurface ice 1136 ! Now we add the flux from the subsurface ice : rho Cd,h U*(1/1+R) * (q1-qsat_ssi(Tice)) 1137 1138 call compute_Tice(nsoil, ptsoil(ig,:,islope), 1139 & ztsrf(ig), 1140 & h2o_ice_depth(ig,islope), 1141 & Tice(ig,islope)) ! compute ice temperature 1108 1142 1109 !EV subsurface ice 1110 IF(h2o_ice_depth(ig,islope) .gt. 0 .and. lag_layer) 1111 & then 1112 zdqsdif_tot(ig)= 1113 & -zqsurf(ig)/subtimestep 1114 zqsurf(ig)=0 1115 1116 DO ik=0,nsoil-2 ! go through all the layers to find the ice locations 1117 IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and. 1118 & (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN 1119 lice = ik+1 1120 EXIT 1121 ENDIF 1122 ENDDO !of subsurface loop 1123 IF (lice .gt. 1) then !calculate the distance from the layers 1124 dist_up(ig,islope)=(h2o_ice_depth(ig,islope) 1125 & -mlayer(lice-1)) 1126 dist_down(ig,islope)=(mlayer(lice) 1127 & -h2o_ice_depth(ig,islope)) 1128 dist_sum(ig,islope)=dist_up(ig,islope) 1129 & +dist_down(ig,islope) 1130 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp 1131 & *ptsoil(ig,lice-1,islope) 1132 & /dist_sum(ig,islope))+ 1133 & (dist_down(ig,islope)*ptsoil(ig,lice,islope) 1134 & /dist_sum(ig,islope)) 1135 ELSE 1136 Tice(ig,islope)=ptsoil(ig,1,islope) 1137 ENDIF 1138 call watersat(1,Tice(ig,1),pplev(ig,1) 1139 & ,qsat2(ig,1)) 1140 qeq(ig,1)=(ztsrf(ig)/Tice(ig,1)) 1141 & *qsat2(ig,1) 1142 resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig,islope) 1143 & /d_coef(ig,islope))) 1144 !write(*,*)'R=',resist(ig,islope) 1145 !write(*,*)'zice=',h2o_ice_depth(ig,islope) 1146 !write(*,*)'D=',d_coef(ig,islope) 1147 !write(*,*)'zcdh=',zcdh(ig) 1148 zb(ig,1)=zb(ig,1)/resist(ig,1) ! change zb to account subsurface ice 1149 !vdifc algorithem !!!!!needs to change reseist io to the mean 1150 !beacuse of the slopes!!!!! 1151 z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) 1152 zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig) 1153 zd(ig,nlay)=zb(ig,nlay)*z1(ig) 1154 DO ilay=nlay-1,2,-1 1155 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ 1156 $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) 1157 zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+ 1158 $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) 1159 zd(ig,ilay)=zb(ig,ilay)*z1(ig) 1160 ENDDO 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)) * z1(ig) 1165 zd(ig,1)=zb(ig,1)*z1(ig) 1166 zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) 1167 zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig,islope) 1168 & *(zq1temp(ig)-qeq(ig,islope)) 1169 call write_output('zdq_ssi', 1170 & '','',zdqsdif_ssi(ig,1)) 1171 call write_output('qeq', 1172 & '','',qeq(ig,1)) 1173 call write_output('q1', 1174 & '','',zq1temp(ig)) 1175 !write(*,*)'qeq=',qeq(ig,1) 1176 !write(*,*)'q1=',zq1temp(ig) 1177 !write(*,*)'zdqsdif_ssi=',zdqsdif_ssi(ig,1) 1178 !I should some all the interactions with the SSI accoridng 1179 !to the statistics and the evaluate 1143 call watersat(1,Tice(ig,islope),pplev(ig,1), 1144 & qsat_ssi(ig,islope)) 1145 1146 qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope) 1147 & *qsat_ssi(ig,islope) 1148 ! Flux from the subsurface 1149 if(old_wsublimation_scheme) then 1150 zdqsdif_ssi_atm(ig,islope)=rho(ig) 1151 & *dryness(ig)*zcdv(ig,islope) 1152 & *1/(1+resist(ig,islope)) 1153 & *(zq1temp(ig)-qsat_ssi(ig,islope)) 1154 else 1155 zdqsdif_ssi_atm(ig,islope)=rho(ig)* 1156 & *dryness(ig) *zcdh(ig,islope) 1157 & *1/(1+resist(ig,islope)) 1158 & *(zq1temp(ig)-qsat_ssi(ig,islope)) 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 1172 1173 else 1174 ! No atm <-> subsurface exchange, we do it the usual way 1175 zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep 1176 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) 1177 zc(ig,1)=(za(ig,1)* 1178 & zq_tmp_vap(ig,1,igcm_h2o_vap)+ 1179 & zb(ig,2)*zc(ig,2) + 1180 & (-zdqsdif_tot(ig)) *subtimestep) *z1(ig) 1181 zq1temp(ig)=zc(ig,1) 1182 endif ! Of subsurface <-> atmosphere exchange 1183 endif ! sublimating more than available frost & surface - frost exchange 1184 endif !if .not.watercaptag(ig) 1185 1186 ! -------------------------------------------------------------------------------------------------------------------------------- 1187 ! We check possible frost subsurface ice interaction: since there is no subsurface water ice mass reservoir represented (yet), 1188 ! we do not include their effect on the mass of surface frost. 1189 ! -------------------------------------------------------------------------------------------------------------------------------- 1190 1191 if((h2o_ice_depth(ig,islope).gt.0).and.lag_layer 1192 & .and.(.not.adsorption_soil)) then 1193 ! First case: still frost at the surface but no watercaptag 1194 if(((watercaptag(ig))).or. 1195 & (((-zdqsdif_tot(ig)*subtimestep) 1196 & .lt.(zqsurf(ig))) 1197 & .and. (zqsurf(ig).gt.tol_frost))) then 1198 ! Still frost at the surface: we consider the possibility to have subsurface <-> frost exchange 1199 ! The flux between frost and ssi is D/zice *(qsat(Tsurf)-qsat_ssi(Tice)) 1200 call compute_Tice(nsoil, ptsoil(ig,:,islope), 1201 & ztsrf(ig), 1202 & h2o_ice_depth(ig,islope), 1203 & Tice(ig,islope)) ! compute ice temperature 1204 1205 call watersat(1,Tice(ig,islope),pplev(ig,1), 1206 & qsat_ssi(ig,islope)) 1207 1208 qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope) 1209 & *qsat_ssi(ig,islope) 1210 1211 zdqsdif_ssi_frost(ig,islope)= d_coef(ig,islope) 1212 & /h2o_ice_depth(ig,islope) 1213 & *rho(ig)*dryness(ig) 1214 & *(qsat(ig)-qsat_ssi(ig,islope)) 1215 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) 1217 zdqsdif_tot(ig) = zdqsdif_tot(ig) - 1218 & zdqsdif_ssi_frost(ig,islope) 1219 endif ! watercaptag or frost at the surface 1220 endif ! interaction frost <-> subsurface ice 1180 1221 1181 !write(*,*) "zdq_ssi*t", zdqsdif_ssi(ig)*subtimestep1182 if (zdqsdif_ssi(ig,1)<0) then1183 !zdqsdif_surf(ig)=zdqsdif_ssi(ig,1)-(zqsurf(ig)/subtimestep)1184 zdqsdif_tot(ig)=zdqsdif_tot(ig)+zdqsdif_ssi(ig,1)1185 call write_output('zdq_zdqssi',1186 & '','',zdqsdif_tot(ig)+zdqsdif_ssi(ig,1))1187 endif1188 ELSE1189 c pdqsdif > 0 : ice condensing1190 c pdqsdif < 0 : ice subliming1191 c write(*,*) "subliming more than available frost: qsurf!"1192 zdqsdif_tot(ig)=1193 & -zqsurf(ig)/subtimestep1194 c write(*,*)'flux vers le sol=',pdqsdif(ig,nq)1195 z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2)))1196 zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,igcm_h2o_vap)+1197 $ zb(ig,2)*zc(ig,2) +1198 $ (-zdqsdif_tot(ig)) *subtimestep) *z1(ig)1199 zq1temp(ig)=zc(ig,1)1200 ENDIF !if h2o_ice_depth>0 and lag_layer1201 endif !if .not.watercaptag(ig)1202 endif ! if sublim more than surface1203 1204 ! subsurface ice <--> frost interaction !EV1205 IF(h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer1206 & .and. zqsurf(ig) .gt. 0) then1207 DO ik=0,nsoil-2 ! go through all the layers to find the ice locations1208 IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.1209 & (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN1210 lice = ik+11211 EXIT1212 ENDIF1213 ENDDO !of subsurface loop1214 IF (lice .gt. 1) then !calculate the distance from the layers1215 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)1216 & -mlayer(lice-1))1217 dist_down(ig,islope)=(mlayer(lice)1218 & -h2o_ice_depth(ig,islope))1219 dist_sum(ig,islope)=dist_up(ig,islope)1220 & +dist_down(ig,islope)1221 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp1222 & *ptsoil(ig,lice-1,islope)1223 & /dist_sum(ig,islope))+1224 & (dist_down(ig,islope)*ptsoil(ig,lice,islope)1225 & /dist_sum(ig,islope))1226 ELSE1227 Tice(ig,islope)=ptsoil(ig,1,islope)1228 ENDIF1229 1230 call watersat(1,Tice(ig,1),pplev(ig,1)1231 & ,qsat2(ig,1))1232 qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))1233 & *qsat2(ig,1)1234 ! write(*,*)'icedep=',h2o_ice_depth(ig,1)1235 ! write(*,*)'qeq=',qeq(ig,1)1236 ! write(*,*)'d=',d_coef(ig,1)1237 ! write(*,*)'qsat=',qsat(ig)1238 ! write(*,*)'dry=',dryness(ig)1239 ! write(*,*)'rho=',rho(ig)1240 zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)1241 & /h2o_ice_depth(ig,1))1242 & *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))1243 !needs to change to the mean of eq1244 zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)1245 !!!!! zdsqdif_ssi_frosst need to be changed to an1246 !average1247 ELSEIF (h2o_ice_depth(ig,islope) .gt. 4e-4 .and. lag_layer1248 & .and. watercaptag(ig)) then1249 DO ik=0,nsoil-2 ! go through all the layers to find the ice locations1250 IF((mlayer(ik).le.h2o_ice_depth(ig,islope)).and.1251 & (mlayer(ik+1).gt.h2o_ice_depth(ig,islope))) THEN1252 lice = ik+11253 EXIT1254 ENDIF1255 ENDDO !of subsurface loop1256 IF (lice .gt. 1) then !calculate the distance from the layers1257 dist_up(ig,islope)=(h2o_ice_depth(ig,islope)1258 & -mlayer(lice-1))1259 dist_down(ig,islope)=(mlayer(lice)1260 & -h2o_ice_depth(ig,islope))1261 dist_sum(ig,islope)=dist_up(ig,islope)1262 & +dist_down(ig,islope)1263 Tice(ig,islope)=(dist_up(ig,islope) ! Linear interp to calculate the temp1264 & *ptsoil(ig,lice-1,islope)1265 & /dist_sum(ig,islope))+1266 & (dist_down(ig,islope)*ptsoil(ig,lice,islope)1267 & /dist_sum(ig,islope))1268 ELSE1269 Tice(ig,islope)=ptsoil(ig,1,islope)1270 ENDIF1271 1272 call watersat(1,Tice(ig,1),pplev(ig,1)1273 & ,qsat2(ig,1))1274 qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))1275 & *qsat2(ig,1)1276 1277 zdqsdif_ssi_frost(ig,1)=(d_coef(ig,1)1278 & /h2o_ice_depth(ig,1))1279 & *rho(ig)*dryness(ig)*(qsat(ig)-qeq(ig,1))1280 zdqsdif_tot(ig)=zdqsdif_tot(ig)-zdqsdif_ssi_frost(ig,1)1281 !needs to change to the mean of eq1282 ENDIF1283 ! call write_output('subtimestep',1284 ! & 'vdifc substimestep length','s',subtimestep)1285 ! ENDDO !subsurface ice subslope1286 1287 1288 1222 1289 1223 c Starting upward calculations for water : 1290 c Actualisation de h2o_vap dans le premier niveau1291 1224 zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig) 1292 1225 c Take into account the H2O latent heat impact on the surface temperature … … 1307 1240 zqsurf(ig)= zqsurf(ig)+( 1308 1241 & zdqsdif_tot(ig))*subtimestep 1309 if (zqsurf(ig)<0 .and. 1310 & (.not.watercaptag(ig))) then 1242 if (zqsurf(ig)<0 .and.(.not.watercaptag(ig))) then 1311 1243 zqsurf(ig)=0 1312 1244 endif 1313 1245 zdqsdif_ssi_atm_tot(ig,islope) = 1246 & zdqsdif_ssi_atm_tot(ig,islope) 1247 & + zdqsdif_ssi_atm(ig,islope) 1248 zdqsdif_ssi_frost_tot(ig,islope) = 1249 & zdqsdif_ssi_frost_tot(ig,islope) 1250 & + zdqsdif_ssi_frost(ig,islope) 1314 1251 c Monitoring instantaneous latent heat flux in W.m-2 : 1315 1252 zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+ … … 1326 1263 endif 1327 1264 1328 c Fin du sous pas de temps1265 c End of the subtimestep 1329 1266 ENDDO ! tsub=1,nsubtimestep 1330 1267 … … 1339 1276 & cos(pi*def_slope_mean(islope)/180.)) 1340 1277 & /ptimestep 1278 1279 zdqsdif_ssi_tot(ig,islope) = 1280 & zdqsdif_ssi_atm_tot(ig,islope) 1281 & + zdqsdif_ssi_frost_tot(ig,islope) 1341 1282 c if subliming more than qsurf(ice) and on watercaptag, water 1342 1283 c sublimates from watercap reservoir (dwatercap_dif is <0) … … 1459 1400 & "Ground ice latent heat flux", 1460 1401 & "W.m-2",surf_h2o_lh(:,iflat)) 1461 ! call write_output('zdqsdif_ssi_frost', 1462 ! & 'Flux between frost and subsurface','kg.m-2.s-1', 1463 ! & zdqsdif_ssi_frost(:,1)) 1464 1465 ! call write_output('zdq_subtimestep', 1466 ! & 'Actual flux zdqsdif_surf*subtimestep', 1467 ! & 'kg.m-2',zdqsdif_tot(:)*subtimestep) 1468 ! call write_output('zdq_end', 1469 ! & 'Flux after all contributions', 1470 ! & 'kg.m-2.s-1',zdqsdif_tot(:)) 1402 1403 call write_output('zdqsdif_ssi_frost_tot', 1404 & 'Flux between frost and subsurface ice','kg.m-2.s-1', 1405 & zdqsdif_ssi_frost_tot(:,iflat)) 1406 1407 call write_output('zdqsdif_ssi_atm_tot', 1408 & 'Flux between atmosphere and subsurface ice','kg.m-2.s-1', 1409 & zdqsdif_ssi_atm_tot(:,iflat)) 1410 1411 call write_output('zdqsdif_ssi_tot', 1412 & 'Total flux echange with subsurface ice','kg.m-2.s-1', 1413 & zdqsdif_ssi_tot(:,iflat)) 1414 1415 1471 1416 C Diagnostic output for HDO 1472 1417 ! if (hdo) then … … 1553 1498 1554 1499 END SUBROUTINE make_tsub 1500 1501 1502 c==================================== 1503 1504 SUBROUTINE compute_Tice(nsoil, ptsoil, ptsurf, ice_depth, Tice) 1505 1506 c Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells. 1507 use comsoil_h, only: layer, mlayer 1508 1509 implicit none 1510 integer,intent(in) :: nsoil ! Number of soil layers 1511 real,intent(in) :: ptsoil(nsoil) ! Soil temperature (K) 1512 real,intent(in) :: ptsurf ! Soil temperature (K) 1513 real,intent(in) :: ice_depth ! Ice depth (m) 1514 real,intent(out) :: Tice ! Ice temperatures (K) 1515 1516 c Local 1517 integer :: ik ! Loop variables 1518 integer :: indexice ! Index of the ice 1519 1520 c Code: 1521 indexice = -1 1522 if(ice_depth.lt.mlayer(0)) then 1523 indexice = 0. 1524 else 1525 do ik = 0,nsoil-2 ! go through all the layers to find the ice locations 1526 if((mlayer(ik).le.ice_depth).and. 1527 & (mlayer(ik+1).gt.ice_depth)) then 1528 indexice = ik+1 1529 exit 1530 endif 1531 enddo 1532 endif 1533 1534 if(indexice.lt.0) then 1535 call abort_physic("vdifc - compute Tice", 1536 & "subsurface ice is below the last soil layer",1) 1537 else 1538 if(indexice .ge. 1) then ! Linear inteprolation between soil temperature 1539 Tice = (ptsoil(indexice)-ptsoil(indexice+1)) 1540 & /(mlayer(indexice-1)-mlayer(indexice)) 1541 & *(ice_depth-mlayer(indexice)) + ptsoil(indexice+1) 1542 else ! Linear inteprolation between the 1st soil temperature and the surface temperature 1543 Tice = (ptsoil(1) - ptsurf)/mlayer(0) 1544 & *(ice_depth-mlayer(0)) + ptsoil(1) 1545 endif ! index ice >=0 1546 endif !indexice <0 1547 1548 1549 END SUBROUTINE compute_Tice 1555 1550 END MODULE vdifc_mod
Note: See TracChangeset
for help on using the changeset viewer.
