Changeset 2885 for trunk/LMDZ.COMMON/libf/evolution/pem.F90
- Timestamp:
- Jan 31, 2023, 9:52:30 PM (22 months ago)
- File:
-
- 1 edited
-
trunk/LMDZ.COMMON/libf/evolution/pem.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r2863 r2885 37 37 albedodat, zmea, zstd, zsig, zgam, zthe, & 38 38 hmons, summit, base,albedo_h2o_frost, & 39 frost_albedo_threshold,emissiv 39 frost_albedo_threshold,emissiv,watercaptag 40 40 use dimradmars_mod, only: totcloudfrac, albedo 41 41 use dust_param_mod, only: tauscaling … … 90 90 TI_PEM,inertiedat_PEM,alph_PEM, beta_PEM, & ! soil thermal inertia 91 91 tsoil_PEM, mlayer_PEM,layer_PEM, & !number of subsurface layers, soil mid layer depths 92 fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys ! geothermal flux, mass of co2 & h2o in the regolith 92 fluxgeo, co2_adsorbded_phys, h2o_adsorbded_phys, & ! geothermal flux, mass of co2 & h2o in the regolith 93 water_reservoir 93 94 use adsorption_mod, only : regolith_adsorption 94 95 … … 187 188 LOGICAL :: STOPPING_co2 ! Logical : is the criterion (% of change in the surface of sublimating water ice) reached? 188 189 LOGICAL :: STOPPING_1_co2 ! Logical : is there still water ice to sublimate? 190 LOGICAL :: STOPPING_pressure 191 INTEGER :: criterion_stop !Which criterion is reached : 1=surf h20, 2=surf_co2, 3=ps, 4=orb.param 192 189 193 190 194 REAL, dimension(:,:,:),allocatable :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg] … … 275 279 INTEGER :: year_iter_max ! maximum number of iterations before stopping 276 280 REAL :: timestep ! timestep [s] 281 REAL WC_sum 282 277 283 #ifdef CPP_STD 278 284 ! INTEGER :: nsplit_phys=1 … … 418 424 enddo 419 425 enddo 426 427 call surfini(ngrid,qsurf) 420 428 421 429 #else … … 793 801 watersurf_density_phys_ave,watersoil_density_phys_PEM_ave, & 794 802 co2_adsorbded_phys,delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded) 803 804 do ig=1,ngrid 805 do islope=1,nslope 806 qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+watercap_slope(ig,islope)+water_reservoir(ig) 807 enddo 808 enddo 795 809 796 810 if(soil_pem) then … … 1139 1153 call criterion_ice_stop_water_slope(cell_area,ini_surf_h2o,qsurf_slope(:,igcm_h2o_ice,:),STOPPING_water,ngrid,initial_h2o_ice_slope) 1140 1154 1141 call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2, ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope)1155 call criterion_ice_stop_slope(cell_area,ini_surf_co2,co2ice_slope,STOPPING_co2,STOPPING_pressure,ngrid,initial_co2_ice_sublim_slope,global_ave_press_GCM,global_ave_press_new,nslope) 1142 1156 1143 1157 year_iter=year_iter+dt_pem … … 1146 1160 if (STOPPING_water) then 1147 1161 print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water 1162 criterion_stop=1 1163 endif 1164 if (STOPPING_1_water) then 1165 print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water 1166 criterion_stop=1 1167 endif 1168 if (STOPPING_co2) then 1169 print *, "STOPPING because surface of co2 ice sublimating is too low, see message above", STOPPING_co2 1170 criterion_stop=2 1171 endif 1172 if (STOPPING_1_co2) then 1173 print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2 1174 criterion_stop=2 1175 endif 1176 if (STOPPING_pressure) then 1177 print *, "STOPPING because surface global pressure changed too much, see message above", STOPPING_pressure 1178 criterion_stop=3 1148 1179 endif 1149 1180 if (year_iter.ge.year_iter_max) then 1150 1181 print *, "STOPPING because maximum number of iterations reached" 1182 criterion_stop=4 1151 1183 endif 1152 if (STOPPING_1_water) then 1153 print *, "STOPPING because tendencies on water ice=0, see message above", STOPPING_1_water 1154 endif 1155 if (STOPPING_co2) then 1156 print *, "STOPPING because surface of co2 ice sublimating is too low or global pressure changed too much, see message above", STOPPING_co2 1157 endif 1158 if (STOPPING_1_co2) then 1159 print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2 1160 endif 1161 1162 1163 1164 if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2) then 1184 1185 1186 1187 if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2 .or. STOPPING_pressure) then 1165 1188 exit 1166 1189 else … … 1197 1220 qsurf(ig,igcm_co2_ice)=co2ice(ig) 1198 1221 #endif 1199 ENDDO ! of DO ig=1,ngrid1200 ! H2o ice1201 DO ig = 1,ngrid1202 qsurf(ig,igcm_h2o_ice) = 0.1203 DO islope = 1,nslope1204 qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) &1205 * subslope_dist(ig,islope) / &1206 cos(pi*def_slope_mean(islope)/180.)1207 ENDDO1208 1222 ENDDO ! of DO ig=1,ngrid 1209 1223 … … 1319 1333 enddo 1320 1334 enddo 1335 1336 DO ig=1,ngrid 1337 if(watercaptag(ig)) then 1338 print *, "INNN" 1339 WC_sum=0. 1340 DO islope=1,nslope 1341 if(qsurf_slope(ig,igcm_h2o_ice,islope).GT. (watercap_slope(ig,islope)+water_reservoir(ig))) then 1342 qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)) 1343 else 1344 watercap_slope(ig,islope)=watercap_slope(ig,islope)+qsurf_slope(ig,igcm_h2o_ice,islope)-(watercap_slope(ig,islope)+water_reservoir(ig)/nslope) 1345 qsurf_slope(ig,igcm_h2o_ice,islope)=0. 1346 endif 1347 WC_sum=WC_sum+watercap_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) 1348 watercap_slope(ig,islope)=0. 1349 enddo 1350 water_reservoir(ig)=water_reservoir(ig)+WC_sum 1351 endif 1352 enddo 1353 1354 ! H2o ice 1355 DO ig = 1,ngrid 1356 qsurf(ig,igcm_h2o_ice) = 0. 1357 DO islope = 1,nslope 1358 qsurf(ig,igcm_h2o_ice) = qsurf(ig,igcm_h2o_ice) + qsurf_slope(ig,igcm_h2o_ice,islope) & 1359 * subslope_dist(ig,islope) / & 1360 cos(pi*def_slope_mean(islope)/180.) 1361 ENDDO 1362 ENDDO ! of DO ig=1,ngrid 1363 1364 DO ig=1,ngrid 1365 ! DO islope=1,nslope 1366 if(qsurf(ig,igcm_h2o_ice).GT.500) then 1367 watercaptag(ig)=.true. 1368 DO islope=1,nslope 1369 qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250 1370 water_reservoir(ig)=water_reservoir(ig)+250 1371 ENDDO 1372 endif 1373 ! enddo 1374 enddo 1375 1376 DO ig=1,ngrid 1377 ! DO islope=1,nslope 1378 if(water_reservoir(ig).LT. 10) then 1379 watercaptag(ig)=.false. 1380 qsurf(ig,igcm_h2o_ice)=qsurf(ig,igcm_h2o_ice)+water_reservoir(ig) 1381 DO islope=1,nslope 1382 qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)+water_reservoir(ig) 1383 ENDDO 1384 endif 1385 ! enddo 1386 enddo 1387 1388 DO ig = 1,ngrid 1389 watercap(ig) = 0. 1390 DO islope = 1,nslope 1391 watercap(ig) = watercap(ig) + watercap_slope(ig,islope) & 1392 * subslope_dist(ig,islope) / & 1393 cos(pi*def_slope_mean(islope)/180.) 1394 ENDDO 1395 ENDDO ! of DO ig=1,ngrid 1396 1321 1397 1322 1398 !------------------------ … … 1386 1462 1387 1463 call pemdem1("restartfi_PEM.nc",year_iter,nsoilmx_PEM,ngrid,nslope , & 1388 tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys) 1464 tsoil_PEM, TI_PEM, ice_depth,co2_adsorbded_phys,h2o_adsorbded_phys,water_reservoir) 1465 1466 call info_run_PEM(year_iter, criterion_stop) 1389 1467 1390 1468 print *, "restartfi_PEM.nc has been written"
Note: See TracChangeset
for help on using the changeset viewer.
