Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4157) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 4170) @@ -52,5 +52,5 @@ use surf_ice, only: evolve_co2ice, evolve_h2oice, balance_h2oice_reservoirs use surf_temp, only: tsurf_PCM, adapt_tsurf2disappearedice -use tendencies, only: compute_tendice, evolve_tend_co2ice, evolve_tend_h2oice +use tendencies, only: compute_tendice, evolve_tend_co2ice, evolve_flux_ssice use tracers, only: adapt_tracers2pressure use utility, only: real2str @@ -72,8 +72,8 @@ logical(k4), dimension(:,:), allocatable :: is_co2ice_flow ! Flag for location of CO2 glacier flow logical(k4), dimension(:,:), allocatable :: is_h2oice_flow ! Flag for location of H2O glacier flow -real(dp), dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg.m-2] -real(dp), dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg.m-2] -real(dp), dimension(:,:,:), allocatable :: minPCM_h2ofrost ! Minimum of H2O frost over the last PCM year [kg.m-2] -real(dp), dimension(:,:,:), allocatable :: minPCM_co2frost ! Minimum of CO2 frost over the last PCM year [kg.m-2] +real(dp), dimension(:,:,:), allocatable :: minPCM_h2operice ! Minimum of H2O perennial ice over the last PCM year [kg/m2] +real(dp), dimension(:,:,:), allocatable :: minPCM_co2perice ! Minimum of CO2 perennial ice over the last PCM year [kg/m2] +real(dp), dimension(:,:,:), allocatable :: minPCM_h2ofrost ! Minimum of H2O frost over the last PCM year [kg/m2] +real(dp), dimension(:,:,:), allocatable :: minPCM_co2frost ! Minimum of CO2 frost over the last PCM year [kg/m2] ! Surface-related: real(dp), dimension(:,:), allocatable :: tsurf_avg_yr1 ! Average surface temperature of the second to last PCM run [K] @@ -141,5 +141,5 @@ call load_xios_data(ps_avg,ps_ts,tsurf_avg,tsurf_avg_yr1,tsoil_avg,tsoil_ts,h2o_surfdensity_avg,h2o_soildensity_avg, & - q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost) + q_h2o_ts,q_co2_ts,minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost,flux_ssice_avg) ! Initiate soil settings and TI @@ -169,7 +169,7 @@ call print_msg('> Computing surface ice tendencies',LVL_NFO) call compute_tendice(minPCM_h2operice,minPCM_h2ofrost,h2o_ice,d_h2oice) -call print_msg('H2O ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO) +call print_msg('H2O ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_h2oice))//' | '//real2str(maxval(d_h2oice)),LVL_NFO) call compute_tendice(minPCM_co2perice,minPCM_co2frost,co2_ice,d_co2ice) -call print_msg('CO2 ice tendencies [kg/m2/yr] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO) +call print_msg('CO2 ice tendencies [kg/m2/y] (min|max): '//real2str(minval(d_co2ice))//' | '//real2str(maxval(d_co2ice)),LVL_NFO) deallocate(minPCM_h2operice,minPCM_co2perice,minPCM_h2ofrost,minPCM_co2frost) @@ -220,6 +220,4 @@ ! Conversion to surface ice call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth) - ! We put the sublimating tendency coming from subsurface ice into the overall tendency - !where (h2oice_depth > 0. .and. zdqsdif_ssi_tot < -minieps) d_h2oice = zdqsdif_ssi_tot end if call allocate_loop_state() @@ -227,5 +225,5 @@ idt = 0 do while (n_yr_run < min(nmax_yr_runorb,nmax_yr_run) .and. n_yr_sim < ntot_yr_sim) - call print_msg('**** Iteration of the PEM run [plnt yr]: '//real2str(n_yr_run + dt),LVL_NFO) + call print_msg('**** Iteration of the PEM run [plnt y]: '//real2str(n_yr_run + dt),LVL_NFO) ! Evolve global surface pressure call evolve_pressure(d_co2ice,delta_co2_ads,do_sorption,ps_avg_glob_old,ps_avg_glob,ps_avg) @@ -237,5 +235,8 @@ allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) if (do_layering) then + allocate(d_h2oice_new(ngrid,nslope)) h2oice_depth_old(:,:) = h2oice_depth(:,:) + ! The sublimating tendency coming from subsurface ice is given through the surface ice tendency + where (h2oice_depth(:,:) > 0. .and. flux_ssice_avg(:,:) < 0._dp) d_h2oice(:,:) = flux_ssice_avg(:,:) do islope = 1,nslope do i = 1,ngrid @@ -247,5 +248,4 @@ call layering2surfice(layerings_map,h2o_ice,co2_ice,h2oice_depth) ! Balance H2O ice reservoirs - allocate(d_h2oice_new(ngrid,nslope)) call stopping_crit_h2o(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,stopcrit) call balance_h2oice_reservoirs(S_atm_2_h2o,S_h2o_2_atm,S_atm_2_h2oice,S_h2oice_2_atm,h2o_ice,d_h2oice,d_h2oice_new) @@ -253,4 +253,5 @@ else zlag(:,:) = 0._dp + zshift_surf(:,:) = 0._dp call evolve_h2oice(delta_h2o_ads,delta_icetable,h2o_ice,d_h2oice,zshift_surf,stopcrit) call evolve_co2ice(co2_ice,d_co2ice,zshift_surf) @@ -306,8 +307,8 @@ write(num,'("_slope",i2.2)') islope end if - call write_diagevo('h2oice'//trim(num),'H2O ice','kg.m-2',h2o_ice(:,islope),(/dim_ngrid/)) - call write_diagevo('co2ice'//trim(num),'CO2 ice','kg.m-2',co2_ice(:,islope),(/dim_ngrid/)) - call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg.m-2.yr-1',d_h2oice(:,islope),(/dim_ngrid/)) - call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg.m-2.yr-1',d_co2ice(:,islope),(/dim_ngrid/)) + call write_diagevo('h2oice'//trim(num),'H2O ice','kg/m2',h2o_ice(:,islope),(/dim_ngrid/)) + call write_diagevo('co2ice'//trim(num),'CO2 ice','kg/m2',co2_ice(:,islope),(/dim_ngrid/)) + call write_diagevo('d_h2oice'//trim(num),'H2O ice tendency','kg/m2/y',d_h2oice(:,islope),(/dim_ngrid/)) + call write_diagevo('d_co2ice'//trim(num),'CO2 ice tendency','kg/m2/y',d_co2ice(:,islope),(/dim_ngrid/)) if (co2ice_flow) then call write_diagevo('Flow_co2ice'//trim(num),'CO2 ice flow location','T/F',merge(1._dp,0._dp,is_co2ice_flow(:,islope)),(/dim_ngrid/)) @@ -328,6 +329,6 @@ call write_diagevo('inertiesoil'//trim(num),'Thermal inertia','SI',TI(:,:,islope),(/dim_ngrid,dim_nsoil/)) if (do_sorption) then - call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg.m-2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) - call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg.m-2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('co2_ads_reg'//trim(num),'CO2 adsorbed in regolith','kg/m2',co2_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) + call write_diagevo('h2o_ads_reg'//trim(num),'H2O adsorbed in regolith','kg/m2',h2o_ads_reg(:,:,islope),(/dim_ngrid,dim_nsoil/)) end if end if @@ -361,17 +362,9 @@ ! Evolve the tendencies call evolve_tend_co2ice(d_co2ice_ini,co2_ice,emissivity_PCM,q_co2_ts_ini,q_co2_ts,ps_ts,ps_avg_glob_ini,ps_avg_glob,d_co2ice) -!~ if (do_layering) then -!~ do i = 1,ngrid -!~ do islope = 1,nslope -!~ if (is_h2oice_sublim_ini(i,islope) .and. h2oice_depth(i,islope) > 0._dp) call evolve_tend_h2oice(h2oice_depth_old(i,islope),h2oice_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope)) -!~ end do -!~ end do -! else -! do i = 1,ngrid -! do islope = 1,nslope -! call evolve_tend_h2oice(icetable_depth_old(i,islope),icetable_depth(i,islope),tsurf_avg(i,islope),tsoil_ts_old(i,:,islope,:),tsoil_ts(i,:,islope,:),d_h2oice(i,islope)) -! end do -! end do -!~ end if + if (do_layering) then + call evolve_flux_ssice(h2oice_depth_old,h2oice_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg) + else if (icetable_equilibrium .or. icetable_dynamic) then + call evolve_flux_ssice(icetable_depth_old,icetable_depth,tsurf_avg,tsoil_ts_old,tsoil_ts,flux_ssice_avg) + end if ! Increment time @@ -383,5 +376,5 @@ if (backup_rate > 0) then if (mod(idt,backup_rate) == 0) call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, & - icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map,idt) + icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map,idt) end if @@ -415,5 +408,5 @@ call save_clim_state(h2o_ice,co2_ice,tsurf_avg,tsurf_dev,tsoil_avg,tsoil_dev,ps_avg,ps_dev,ps_avg_glob,ps_avg_glob_ini, & - icetable_depth,icetable_thickness,ice_porefilling,h2o_ads_reg,co2_ads_reg,layerings_map) + icetable_depth,icetable_thickness,ice_porefilling,h2oice_depth,flux_ssice_avg,h2o_ads_reg,co2_ads_reg,layerings_map) ! Update the duration information of the workflow