Index: trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90 (revision 2849) @@ -32,5 +32,5 @@ REAL :: rho_regolith = 2000. ! density of the reoglith, Buhler & Piqueux 2021 real :: alpha_clapeyron = -6143.7! eq. (2) in Murphy & Koop 2005 - real :: beta_clapeyron = 29.9074 ! eq. (2) in Murphy & Koop 2005 + real :: beta_clapeyron = 28.9074 ! eq. (2) in Murphy & Koop 2005 real :: mi = 2.84e-7 ! Mass of h2o per m^2 absorbed Jackosky et al. 1997 real :: as = 18.9e3 ! Specific area, Buhler & Piqueux 2021 @@ -88,6 +88,7 @@ do islope = 1,nslope do iloop = 1,n_1km + K = Ko*exp(e/tsoil_PEM(ig,iloop,islope)) if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then - K = Ko*exp(e/tsoil_PEM(ig,iloop,islope)) + theta_h2o_adsorbded(ig,iloop,islope) = (K*pvapor_avg(ig)/(1+K*pvapor_avg(ig)))**mu m_h2o_adsorbed(ig,iloop,islope) = as*theta_h2o_adsorbded(ig,iloop,islope)*mi*rho_regolith Index: trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/computeice_table.F90 (revision 2849) @@ -1,3 +1,3 @@ - SUBROUTINE computeice_table(timelen,ngrid,nslope,nsoil_GCM,nsoil_PEM,tsoil,tsurf,q_co2,q_h2o,ps,ice_table) + SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table) #ifndef CPP_STD USE comsoil_h, only: inertiedat, volcapa @@ -6,127 +6,29 @@ implicit none - integer,intent(in) :: timelen,ngrid,nslope,nsoil_PEM,nsoil_GCM - real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope,timelen) ! soil temperature, timeseries [K] - real,intent(in) :: tsurf(ngrid,nslope,timelen) ! surface temperature [K] - real,intent(in) :: q_co2(ngrid,timelen) ! MMR tracer co2 [kg/kg] - real,intent(in) :: q_h2o(ngrid,timelen) ! MMR tracer h2o [kg/kg] - real,intent(in) :: ps(ngrid,timelen) ! surface pressure [Pa] - real,intent(out) :: ice_table(ngrid,nslope) ! ice table [m] + integer,intent(in) :: ngrid,nslope,nsoil_PEM + real,intent(in) :: rhowatersurf_ave(ngrid,nslope) ! surface temperature, timeseries [K] + real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope) ! soil water density, yearly average [kg/m^3] + real,intent(inout) :: ice_table(ngrid,nslope) ! ice table [m] + real :: z1,z2 + integer ig, islope,isoil + real :: diff_rho(nsoil_PEM) ! difference of vapor content - real :: m_h2o = 18.01528E-3 - real :: m_co2 = 44.01E-3 - real :: m_noco2 = 33.37E-3 - real :: A,B,z1,z2 - real :: alpha = -6143.7 - real :: beta = 29.9074 - - integer ig, islope,isoil,it - real,allocatable :: mass_mean(:,:) ! mean mass above the surface - real,allocatable :: zplev_mean(:,:) ! pressure above the surface - real,allocatable :: pvapor(:,:) ! partial pressure above the surface - - real,allocatable :: rhovapor(:,:,:) - real,allocatable :: rhovapor_avg(:,:) ! mean vapor_density at the surface yearly averaged - - real :: psv_surf - real,allocatable :: rho_soil(:,:,:,:) ! water vapor in the soil - real,allocatable :: rho_soil_avg(:,:,:) ! water vapor in the soil yearly averaged - - real,allocatable :: diff_rho(:,:,:) ! difference of vapor content - - allocate(rhovapor(ngrid,nslope,timelen)) - allocate(rhovapor_avg(ngrid,nslope)) - allocate(pvapor(ngrid,timelen)) - - allocate(mass_mean(ngrid,timelen)) - allocate(zplev_mean(ngrid,timelen)) - -! 0. Some initializations - - A =(1/m_co2 - 1/m_noco2) - B=1/m_noco2 -! 1. Compute rho surface yearly averaged - -! 1.1 Compute the partial pressure of vapor -!a. the molecular mass into the column - do ig = 1,ngrid - mass_mean(ig,:) = 1/(A*q_co2(ig,:) +B) - enddo - -! b. pressure level - do it = 1,timelen - do ig = 1,ngrid - zplev_mean(ig,it) = ap(1) + bp(1)*ps(ig,it) - enddo - enddo - -! c. Vapor pressure - pvapor(:,:) = mass_mean(:,:)/m_h2o*q_h2o(:,:)*zplev_mean(:,:) - - deallocate(mass_mean) - deallocate(zplev_mean) - -! 1.2 Check if there is frost at the surface and then compute the density -! at the surface - do ig = 1,ngrid - do islope = 1,nslope - do it = 1,timelen - psv_surf = exp(alpha/tsurf(ig,islope,it) +beta) - if ((isnan(psv_surf)).or.(isnan(pvapor(ig,it)))) then - write(*,*) 'pb vapor',ig,islope,it - stop - endif - rhovapor(ig,islope,it) = min(psv_surf,pvapor(ig,it))/tsurf(ig,islope,it) - enddo - enddo - enddo - deallocate(pvapor) - -! 1.3 Density at the surface yearly averaged - rhovapor_avg(:,:) = SUM(rhovapor(:,:,:),3)/timelen - - deallocate(rhovapor) - -! 2. Compute rho soil vapor - - allocate(rho_soil_avg(ngrid,nslope,nsoil_PEM)) - allocate(rho_soil(ngrid,nslope,nsoil_PEM,timelen)) do ig = 1,ngrid - do islope = 1,nslope - do isoil = 1,nsoil_PEM - do it = 1,timelen - rho_soil(ig,islope,isoil,it) = exp(alpha/tsoil(ig,isoil,islope,it) +beta)/tsoil(ig,isoil,islope,it) - enddo - enddo - enddo - enddo - - rho_soil_avg(:,:,:) = SUM( rho_soil(:,:,:,:),4)/timelen - deallocate(rho_soil) - -! 3. Computing ice table - - ice_table (:,:) = -1e4 - - allocate(diff_rho(ngrid,nslope,nsoil_PEM)) + do islope = 1,nslope + ice_table(ig,islope) = -10. do isoil = 1,nsoil_PEM - diff_rho(:,:,isoil) = rhovapor_avg(:,:) - rho_soil_avg(:,:,isoil) -! write(*,*) 'diff =',ig,islope,isoil,diff_rho(ig,islope,isoil),rhovapor_avg(ig,islope) ,rho_soil_avg(ig,islope,isoil) + diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope) + enddo - - deallocate(rhovapor_avg) - deallocate(rho_soil_avg) - do ig = 1,ngrid - do islope = 1,nslope - if(diff_rho(ig,islope,1) > 0) then + if(diff_rho(1) > 0) then ice_table(ig,islope) = 0. else do isoil = 1,nsoil_PEM -1 - if((diff_rho(ig,islope,isoil).lt.0).and.(diff_rho(ig,islope,isoil+1).gt.0.)) then - z1 = (diff_rho(ig,islope,isoil) - diff_rho(ig,islope,isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1)) - z2 = -layer_PEM(isoil+1)*z1 + diff_rho(ig,islope,isoil+1) + if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then + z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1)) + z2 = -layer_PEM(isoil+1)*z1 + diff_rho(isoil+1) ice_table(ig,islope) = -z2/z1 exit @@ -136,6 +38,6 @@ enddo enddo + - deallocate(diff_rho) !======================================================================= Index: trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 (revision 2849) @@ -37,4 +37,6 @@ !======================================================================= + + STOPPING=.false. pos_tend=0. Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 2849) @@ -78,5 +78,5 @@ co2iceflow, beta_slope, capcal_slope,& albedo_slope,emiss_slope,qsurf_slope,& - iflat,ini_comslope_h + iflat,major_slope,ini_comslope_h use time_phylmdz_mod, only: daysec,dtphys USE comconst_mod, ONLY: rad,g,r,cpp,pi @@ -272,8 +272,15 @@ REAL,ALLOCATABLE :: alph_locslope(:,:) REAL,ALLOCATABLE :: beta_locslope(:,:) - + REAL,ALLOCATABLE :: watersurf_density_timeseries(:,:,:,:) + REAL,ALLOCATABLE :: watersoil_density_timeseries(:,:,:,:,:) + REAL,ALLOCATABLE :: watersurf_density_phys_timeseries(:,:,:) + REAL,ALLOCATABLE :: watersurf_density_phys_ave(:,:) + REAL,ALLOCATABLE :: watersoil_density_phys_PEM_timeseries(:,:,:,:) + REAL,ALLOCATABLE :: watersoil_density_phys_PEM_ave(:,:,:) REAL,ALLOCATABLE :: Tsurfave_before_saved(:,:) REAL, ALLOCATABLE :: delta_co2_adsorbded(:) REAL :: totmass_adsorbded + real :: alpha_clap_h2o = -6143.7 + real :: beta_clap_h2o = 28.9074 #ifdef CPP_STD @@ -289,5 +296,5 @@ ! Loop variable LOGICAL :: bool_sublim - INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop + INTEGER :: i,j,ig0,l,ig,nnq,t,islope,ig_loop,islope_loop,iloop,isoil REAL :: beta = 3182.48 REAL :: alpha = 23.3494 @@ -406,5 +413,5 @@ watercap,inertiesoil,nslope,tsurf_slope, & tsoil_slope,co2ice_slope,def_slope,def_slope_mean, & - subslope_dist,albedo_slope,emiss_slope, TI_GCM_start, & + subslope_dist,major_slope,albedo_slope,emiss_slope, TI_GCM_start, & qsurf_slope,watercap_slope) @@ -551,9 +558,15 @@ allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) allocate(delta_co2_adsorbded(ngrid)) - + allocate(watersurf_density_timeseries(iim+1,jjm+1,nslope,timelen)) + allocate(watersoil_density_timeseries(iim+1,jjm+1,nsoilmx,nslope,timelen)) + allocate(watersurf_density_phys_timeseries(ngrid,nslope,timelen)) + allocate(watersurf_density_phys_ave(ngrid,nslope)) + allocate(watersoil_density_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) + allocate(watersoil_density_phys_PEM_ave(ngrid,nsoilmx_PEM,nslope)) print *, "Downloading data Y1..." call read_data_GCM("data_GCM_Y1.nc", min_h2o_ice_s_1,min_co2_ice_s_1,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM_yr1,timelen,min_co2_ice_slope_1,min_h2o_ice_slope_1,& - nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope) + nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, & + watersurf_density_timeseries,watersoil_density_timeseries) ! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value @@ -569,5 +582,6 @@ call read_data_GCM("data_GCM_Y2.nc", min_h2o_ice_s_2,min_co2_ice_s_2,iim,jjm,nlayer,vmr_co2_gcm,ps_GCM,timelen,min_co2_ice_slope_2,min_h2o_ice_slope_2, & - nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope) + nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, & + watersurf_density_timeseries,watersoil_density_timeseries) print *, "Downloading data Y2 done" @@ -622,4 +636,5 @@ allocate(ice_depth(ngrid,nslope)) + ice_depth(:,:) = 0. allocate(TI_GCM_phys(ngrid,nsoilmx,nslope)) @@ -636,11 +651,13 @@ deallocate(co2_ice_GCM_slope) deallocate(TI_GCM) - - if(soil_pem) then - deallocate(tsurf_GCM_timeseries) + + + if(soil_pem) then call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,TI_GCM_phys,TI_PEM) - DO islope = 1,nslope + DO t=1,timelen + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersurf_density_timeseries(:,:,islope,t),watersurf_density_phys_timeseries(:,islope,t)) + ENDDO DO l=1,nsoilmx CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,l,islope),tsoil_ave_phys_yr1(:,l,islope)) @@ -648,12 +665,21 @@ DO t=1,timelen CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_GCM_timeseries(:,:,l,islope,t),tsoil_phys_PEM_timeseries(:,l,islope,t)) + CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t)) ENDDO + ENDDO DO l=nsoilmx+1,nsoilmx_PEM CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave_yr1(:,:,nsoilmx,islope),tsoil_ave_phys_yr1(:,l,islope)) CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,tsoil_ave(:,:,nsoilmx,islope),tsoil_PEM(:,l,islope)) + DO t=1,timelen + watersoil_density_phys_PEM_timeseries(:,l,islope,t) = watersoil_density_phys_PEM_timeseries(:,nsoilmx,islope,t) + ENDDO ENDDO ENDDO - + watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen + watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen + deallocate(watersurf_density_timeseries) + deallocate(watersurf_density_phys_timeseries) + deallocate(watersoil_density_timeseries) deallocate(tsoil_ave_yr1) deallocate(tsoil_ave) @@ -789,7 +815,8 @@ !---------------------------- Read the PEMstart --------------------- - call pemetat0(.false.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, & + call pemetat0(.true.,"startfi_PEM.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_ave_phys_yr1,tsoil_PEM,ice_depth, & tsurf_ave_phys_yr1, tsurf_ave_phys, q_co2_PEM_phys, q_h2o_PEM_phys,ps_phys_timeseries,tsurf_phys_GCM_timeseries,tsoil_phys_PEM_timeseries,& - tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:), co2_adsorbded_phys,delta_co2_adsorbded) + tendencies_h2o_ice_phys_slope,tendencies_co2_ice_phys_slope,co2ice_slope,qsurf_slope(:,igcm_h2o_ice,:),global_ave_press_GCM, co2_adsorbded_phys,delta_co2_adsorbded,& + watersurf_density_phys_ave,watersoil_density_phys_PEM_ave) if(soil_pem) then @@ -848,4 +875,5 @@ year_iter=0 print *, "Max_iter_pem", Max_iter_pem + print *, 'year_iter_max',year_iter_max do while (year_iter.LT.year_iter_max) @@ -856,8 +884,15 @@ global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface enddo + enddo + print *, 'Global average pressure old time step',global_ave_press_old + print *, 'Global average pressure new time step',global_ave_press_new + + do i=1,ngrid if(soil_pem) then global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface endif enddo + print *, 'Global average pressure old time step',global_ave_press_old + print *, 'Global average pressure new time step',global_ave_press_new ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption) @@ -974,5 +1009,5 @@ flag_co2flow_mesh(ig) = 1. ENDIF ! co2ice > hmax - ENDIF ! iflattsoil_lope + ENDIF ! iflat ENDDO !islope ENDDO !ig @@ -1032,8 +1067,8 @@ emiss_slope(ig,islope) = emissiv endif - elseif( co2ice_slope(ig,islope).GT. 1E-5) THEN !Put tsurf as tcond co2 + elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-5) )THEN !Put tsurf as tcond co2 ave=0. do t=1,timelen - if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-5) then + if(co2_ice_GCM_phys_slope(ig,islope,t).gt.1e-3) then ave = ave + beta/(alpha-log(vmr_co2_pem_phys(ig,t)*ps_phys_timeseries(ig,t)/100.)) else @@ -1056,4 +1091,5 @@ enddo + if(soil_pem) then @@ -1071,39 +1107,43 @@ do islope = 1,nslope TI_locslope(:,:) = TI_PEM(:,:,islope) - Tsoil_locslope(:,:) = tsoil_PEM(:,:,islope) - Tsurf_locslope(:) = tsurf_ave_phys(:,islope) alph_locslope(:,:) = alph_PEM(:,:,islope) beta_locslope(:,:) = beta_PEM(:,:,islope) - call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) - call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) do t = 1,timelen - tsoil_phys_PEM_timeseries(:,:,islope,t) = tsoil_phys_PEM_timeseries(:,:,islope,t) + (Tsoil_locslope(:,:)- tsoil_PEM(:,:,islope)) - enddo - TI_PEM(:,:,islope) = TI_locslope(:,:) - tsoil_PEM(:,:,islope) = Tsoil_locslope(:,:) - tsurf_ave_phys(:,islope) = Tsurf_locslope(:) - alph_PEM(:,:,islope) = alph_locslope(:,:) - beta_PEM(:,:,islope) = beta_locslope(:,:) - enddo - - print *, "Update of soil temperature done" - -! Check Nan in the time series - do ig = 1,ngrid - do islope = 1,nslope - do iloop = 1,nsoilmx_PEM - if(isnan(tsoil_PEM(ig,iloop,islope))) then - write(*,*) "Problem : There is nan tsoil", ig, iloop, islope, tsoil_PEM(ig,iloop,islope) - stop - endif + Tsoil_locslope(:,:) = tsoil_phys_PEM_timeseries(:,:,islope,t) + Tsurf_locslope(:) = tsurf_phys_GCM_timeseries(:,islope,t) + + call soil_pem_routine(ngrid,nsoilmx_PEM,.true.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) + call soil_pem_routine(ngrid,nsoilmx_PEM,.false.,TI_locslope,timestep/timelen,Tsurf_locslope,Tsoil_locslope,alph_locslope,beta_locslope) + + + tsoil_phys_PEM_timeseries(:,:,islope,t) = Tsoil_locslope(:,:) + do ig = 1,ngrid + do isoil = 1,nsoilmx_PEM + watersoil_density_phys_PEM_timeseries(ig,l,islope,t) = exp(alpha_clap_h2o/Tsoil_locslope(ig,isoil) + beta_clap_h2o)/Tsoil_locslope(ig,isoil) + if(isnan(Tsoil_locslope(ig,isoil))) then + write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil) + endif enddo enddo enddo + alph_PEM(:,:,islope) = alph_locslope(:,:) + beta_PEM(:,:,islope) = beta_locslope(:,:) + enddo + + + tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen + watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen + print *, "Update of soil temperature done" + + deallocate(TI_locslope) + deallocate(Tsoil_locslope) + deallocate(Tsurf_locslope) + deallocate(alph_locslope) + deallocate(beta_locslope) write(*,*) "Compute ice table" ! II_d.3 Update the ice table - call computeice_table(timelen,ngrid,nslope,nsoilmx,nsoilmx_PEM, tsoil_phys_PEM_timeseries,tsurf_phys_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & - ps_phys_timeseries, ice_depth) + call computeice_table(ngrid,nslope,nsoilmx_PEM,watersurf_density_phys_ave,watersoil_density_phys_PEM_ave,ice_depth) print *, "Update soil propreties" @@ -1153,4 +1193,7 @@ if (STOPPING_water) then print *, "STOPPING because surface of water ice sublimating is too low, see message above", STOPPING_water + endif + if (year_iter.ge.year_iter_max) then + print *, "STOPPING because maximum number of iterations reached" endif if (STOPPING_1_water) then @@ -1210,5 +1253,5 @@ ENDDO ! of DO ig=1,ngrid -! III_a.2 Tsurf and Tsoil update (for startfi) +! III_a.2 Tsoil update (for startfi) if(soil_pem) then @@ -1219,11 +1262,4 @@ endif !soil_pem - DO ig = 1,ngrid - tsurf(ig) = 0. - DO islope = 1,nslope - tsurf(ig) = tsurf(ig) + tsurf_slope(ig,islope) & - * subslope_dist(ig,islope) - ENDDO - ENDDO ! of DO ig=1,ngrid #ifndef CPP_STD @@ -1260,4 +1296,15 @@ ENDDO ENDDO + + + DO ig = 1,ngrid + tsurf(ig) = 0. + DO islope = 1,nslope + tsurf(ig) = tsurf(ig) + (emiss_slope(ig,islope)*tsurf_slope(ig,islope))**4 & + * subslope_dist(ig,islope) + ENDDO + tsurf(ig) = tsurf(ig)**(0.25)/emis(ig) + ENDDO ! of DO ig=1,ngrid + #endif Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 2849) @@ -1,9 +1,9 @@ SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, & - tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, m_co2_regolith_phys,deltam_co2_regolith_phys) + tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, & + global_ave_pressure, m_co2_regolith_phys,deltam_co2_regolith_phys, watersurf_ave,watersoil_ave) use iostart_PEM, only: open_startphy, close_startphy, get_field, get_var use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM use comsoil_h, only: volcapa,inertiedat - use ini_soil_mod, only: ini_icetable use soil_evolution_mod, only: soil_pem,soil_pem_CN use adsorption_mod, only : regolith_co2adsorption @@ -32,5 +32,7 @@ real,intent(in) :: waterice(ngrid,nslope) real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope) - + real, intent(in) :: watersurf_ave(ngrid,nslope) ! surface water ice density, yearly averaged (kg/m^3) + real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope) ! surface water ice density, yearly averaged (kg/m^3) + real, intent(in) :: global_ave_pressure ! mean average pressure on the planet ! outputs @@ -59,4 +61,6 @@ real :: co2_ads_prev(ngrid) real :: year_PEM_read + real :: alpha_clap_h2o = -6143.7 + real :: beta_clap_h2o = 28.9074 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -221,4 +225,10 @@ enddo enddo + do isoil = nsoil_GCM+1,nsoil_PEM + do ig = 1,ngrid + watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope) + enddo + enddo + ENDDO @@ -230,16 +240,9 @@ !3. Ice Table - call get_field("ice_table",ice_table,found) - if(.not.found) then - write(*,*)'PEM settings: failed loading ' - write(*,*)'will reconstruct the values of ice table' - - call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope)) - - else - ! update ice table - call computeice_table(timelen,ngrid,nslope,nsoil_PEM,tsoil_inst,tsurf_inst,q_co2,q_h2o, ps_inst, ice_table) - - endif + + + call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table) + call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM) + print *,'PEMETAT0: ICE TABLE DONE' @@ -259,17 +262,11 @@ deltam_co2_regolith_phys(:) = 0. exit + else + call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys) + endif ENDDO -if (found) then - DO iloop = 1,nsoil_GCM - tsoil_tmp_yr1(:,iloop,:) = tsoil_PEM_yr1(:,iloop,:) - - ENDDO - - call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys) - - -endif + print *,'PEMETAT0: CO2 adsorption done ' @@ -364,6 +361,4 @@ !b) Soil temperature do islope = 1,nslope - - write(*,*) "islope=",islope do ig = 1,ngrid kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa @@ -381,5 +376,4 @@ tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km)) enddo -! write(*,*) "ig, islope, T=", ig,islope,tsoil_PEM(ig,:,islope) enddo @@ -389,12 +383,17 @@ enddo enddo + do isoil = nsoil_GCM+1,nsoil_PEM + do ig = 1,ngrid + watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope) + enddo + enddo enddo print *,'PEMETAT0: TSOIL DONE ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !c) Ice table - do islope = 1,nslope - call ini_icetable(timelen,ngrid,nsoil_PEM,TI_PEM, timestep,tsurf_ave(:,islope),tsoil_PEM(:,:,islope),tsurf_inst(:,islope,:), & - tsoil_inst(:,:,islope,:),q_co2,q_h2o,ps_inst,ice_table(:,islope)) - enddo + call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table) + call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90 (revision 2849) @@ -3,5 +3,6 @@ ! SUBROUTINE read_data_GCM(fichnom,min_h2o_ice_s,min_co2_ice_s,iim_input,jjm_input,nlayer,vmr_co2_gcm,ps_GCM,timelen, & - min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope) + min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, & + watersurf_density,watersoil_density) use netcdf, only: nf90_open,NF90_NOWRITE,nf90_noerr,nf90_strerror, & @@ -39,5 +40,4 @@ REAL, INTENT(OUT) :: q_h2o_GCM(iim_input+1,jjm_input+1,timelen) REAL, INTENT(OUT) :: q_co2_GCM(iim_input+1,jjm_input+1,timelen) - REAL, ALLOCATABLE :: q1_co2_GCM(:,:,:) REAL, INTENT(OUT) :: ps_GCM(iim_input+1,jjm_input+1,timelen) @@ -48,4 +48,6 @@ REAL ,INTENT(OUT) :: tsurf_gcm(iim_input+1,jjm_input+1,nslope,timelen) ! Surface temperature of the concatenated file REAL , INTENT(OUT) :: tsoil_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file + REAL , INTENT(OUT) :: watersurf_density(iim_input+1,jjm_input+1,nslope,timelen) ! Soil temperature of the concatenated file + REAL , INTENT(OUT) :: watersoil_density(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Soil temperature of the concatenated file REAL :: TI_gcm(iim_input+1,jjm_input+1,nsoilmx,nslope,timelen) ! Thermal Inertia of the concatenated file @@ -79,5 +81,4 @@ allocate(co2_ice_s(iim+1,jjm+1,timelen)) - allocate(q1_co2_GCM(iim+1,jjm+1,timelen)) allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen)) allocate(h2o_ice_s(iim+1,jjm+1,timelen)) @@ -161,4 +162,22 @@ print *, "Downloading data for inertiesoil_slope done" + + print *, "Downloading data for watersoil_density ..." + +DO islope=1,nslope + write(num,fmt='(i2.2)') islope + call get_var4("Waterdensity_soil_slope"//num,watersoil_density(:,:,:,islope,:)) +ENDDO + + print *, "Downloading data for watersoil_density done" + + print *, "Downloading data for watersurf_density ..." + +DO islope=1,nslope + write(num,fmt='(i2.2)') islope + call get_var3("Waterdensity_surface"//num,watersurf_density(:,:,islope,:)) +ENDDO + + print *, "Downloading data for watersurf_density done" endif !soil_pem @@ -234,4 +253,10 @@ ENDDO + + deallocate(co2_ice_s) + deallocate(h2o_ice_s_slope) + deallocate(h2o_ice_s) + + CONTAINS Index: trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90 (revision 2849) @@ -8,10 +8,13 @@ USE temps_mod_evol, ONLY: year_bp_ini, year_PEM -#ifndef CPP_STD +#ifndef CPP_STD + USE comconst_mod, ONLY: pi USE planete_h, ONLY: e_elips, obliquit, timeperi #else use planete_mod, only: e_elips, obliquit, timeperi + USE comcstfi_mod, only: pi + #endif - USE comcstfi_mod, only: pi + USE lask_param_mod, only: yearlask,oblask,exlask,lsplask, & end_lask_param_mod,last_ilask Index: trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 (revision 2842) +++ trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 (revision 2849) @@ -26,8 +26,6 @@ REAL :: L = 51058. REAL :: inertie_thresold = 800. ! look for ice - REAL :: inertie_co2glaciers = 2120 ! Mellon et al. 2000 REAL :: inertie_averaged = 250 ! Mellon et al. 2000 - REAL :: ice_inertia = 2120 ! Inertia of ice - REAL :: surfaceice_inertia = 800 ! Inertia of ice + REAL :: ice_inertia = 1200 ! Inertia of ice REAL :: P610 = 610.0 REAL :: m_h2o = 18.01528E-3 @@ -47,29 +45,27 @@ ! 0. Initialisation - do ig = 1,ngrid - do islope = 1,nslope - if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then - ispermanent_h2oglaciers(ig,islope) = 1 - else - ispermanent_h2oglaciers(ig,islope) = 0 - endif +! do ig = 1,ngrid +! do islope = 1,nslope +! if((abs(tend_h2oglaciers(ig,islope)).lt.1e-5).and.(abs(waterice(ig,islope)).gt.(1.e-4))) then +! ispermanent_h2oglaciers(ig,islope) = 1 +! else +! ispermanent_h2oglaciers(ig,islope) = 0 +! endif +! +! if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.(1.e-4))) then +! ispermanent_co2glaciers(ig,islope) = 1 +! else +! ispermanent_co2glaciers(ig,islope) = 0 +! endif +! enddo +! enddo - if((abs(tend_co2glaciers(ig,islope)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then - ispermanent_co2glaciers(ig,islope) = 1 - else - ispermanent_co2glaciers(ig,islope) = 0 - endif - enddo - enddo - - ispermanent_co2glaciers(:,:) = 0 - ispermanent_h2oglaciers(:,:) = 0 ! 1.Ice TI feedback - do islope = 1,nslope - call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope), TI_PEM(:,:,islope)) - enddo +! do islope = 1,nslope +! call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope), TI_PEM(:,:,islope)) +! enddo ! 2. Modification of the regolith thermal inertia.