MODULE ini_soil_mod IMPLICIT NONE CONTAINS subroutine ini_icetable(timelen,ngrid,nsoil_PEM, & therm_i, timestep,tsurf_ave,tsoil_ave,tsurf_inst, tsoil_inst,q_co2,q_h2o,ps,ice_table) use vertical_layers_mod, only: ap,bp use comsoil_h_PEM, only: fluxgeo,layer_PEM,inertiedat_PEM use comsoil_h,only: volcapa, nsoilmx implicit none !----------------------------------------------------------------------- ! Author: LL ! Purpose: Compute soil temperature using an implict 1st order scheme ! ! Note: depths of layers and mid-layers, soil thermal inertia and ! heat capacity are commons in comsoil_PEM.h ! A convergence loop is added until equilibrium !----------------------------------------------------------------------- #include "dimensions.h" !#include "dimphys.h" !#include"comsoil.h" !----------------------------------------------------------------------- ! arguments ! --------- ! inputs: integer,intent(in) :: timelen ! Time length in for time-series data integer,intent(in) :: ngrid ! number of (horizontal) grid-points integer,intent(in) :: nsoil_PEM ! number of soil layers in the PEM real,intent(in) :: timestep ! time step real,intent(in) :: tsurf_ave(ngrid) ! surface temperature 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(in) :: tsurf_inst(ngrid,timelen) ! soil (mid-layer) temperature ! outputs: real,intent(inout) :: tsoil_ave(ngrid,nsoil_PEM) ! soil (mid-layer) temperature. real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,timelen) ! soil (mid-layer) temperature real,intent(out) :: ice_table(ngrid) ! ice table [m] real,intent(inout) :: therm_i(ngrid,nsoil_PEM) ! thermal inertia ! local variables: integer ig,isoil,it,k,iref REAL :: error_depth = 2. REAL :: tsoil_saved(nsoil_PEM) integer :: countmax = 20 integer :: countloop REAL :: tol_error = 0.1 REAL :: ice_inertia = 2120. REAL :: alph_PEM(nsoil_PEM-1) REAL :: beta_PEM(nsoil_PEM-1) real :: rhoc(nsoil_PEM) real :: volcapa_ice = 1.43e7 real :: k_soil(nsoil_PEM) real :: d1(nsoil_PEM),d2(nsoil_PEM),d3(nsoil_PEM),re(nsoil_PEM) real :: Tcol_saved(nsoil_PEM) real :: tsoil_prev(nsoil_PEM) real :: tsurf_prev real :: icedepth_prev 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 real,allocatable :: mass_mean(:) ! mean mass above the surface real,allocatable :: zplev(:) ! pressure above the surface real,allocatable :: pvapor(:) ! partial pressure above the surface real,allocatable :: rhovapor(:) real :: 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 A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 ice_table(:) = 1. do ig = 1,ngrid error_depth = 1. countloop = 0 do while(( error_depth.gt.tol_error).and.(countloop.lt.countmax).and.(ice_table(ig).gt.-1e-20)) countloop = countloop +1 Tcol_saved(:) = tsoil_ave(ig,:) !2. Compute ice table ! 2.1 Compute water density at the surface, yearly averaged allocate(mass_mean(timelen)) ! 1.1 Compute the partial pressure of vapor ! a. the molecular mass into the column mass_mean(:) = 1/(A*q_co2(ig,:) +B) ! b. pressure level allocate(zplev(timelen)) do it = 1,timelen zplev(it) = ap(1) + bp(1)*ps(ig,it) enddo ! c. Vapor pressure allocate(pvapor(timelen)) pvapor(:) = mass_mean(:)/m_h2o*q_h2o(ig,:)*zplev(:) deallocate(zplev) deallocate(mass_mean) ! d! Check if there is frost at the surface and then compute the density ! at the surface allocate(rhovapor(timelen)) do it = 1,timelen psv_surf = exp(alpha/tsurf_inst(ig,it) +beta) rhovapor(it) = min(psv_surf,pvapor(it))/tsurf_inst(ig,it) enddo deallocate(pvapor) rhovapor_avg = SUM(rhovapor(:),1)/timelen deallocate(rhovapor) ! 2.2 Compute water density at the soil layer, yearly averaged allocate(rho_soil(timelen)) allocate(rho_soil_avg(nsoil_PEM)) do isoil = 1,nsoil_PEM do it = 1,timelen rho_soil(it) = exp(alpha/tsoil_inst(ig,isoil,it) +beta)/tsoil_inst(ig,isoil,it) enddo rho_soil_avg(isoil) = SUM(rho_soil(:),1)/timelen enddo deallocate(rho_soil) !2.3 Final: compute ice table icedepth_prev = ice_table(ig) ice_table(ig) = -1 allocate(diff_rho(nsoil_PEM)) do isoil = 1,nsoil_PEM diff_rho(isoil) = rhovapor_avg - rho_soil_avg(isoil) enddo deallocate(rho_soil_avg) if(diff_rho(1) > 0) then ice_table(ig) = 0. else do isoil = 1,nsoil_PEM -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) = -z2/z1 exit endif enddo endif deallocate(diff_rho) !3. Update Soil Thermal Inertia if (ice_table(ig).gt. 0.) then if (ice_table(ig).lt. 1e-10) then do isoil = 1,nsoil_PEM therm_i(ig,isoil)=ice_inertia enddo else ! 4.1 find the index of the mixed layer iref=0 ! initialize iref do k=1,nsoil_PEM ! loop on layers if (ice_table(ig).ge.layer_PEM(k)) then iref=k ! pure regolith layer up to here else ! correct iref was obtained in previous cycle exit endif enddo ! 4.2 Build the new ti do isoil=1,iref therm_i(ig,isoil) =inertiedat_PEM(ig,isoil) enddo if (iref.lt.nsoil_PEM) then if (iref.ne.0) then ! mixed layer therm_i(ig,iref+1)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ & (((ice_table(ig)-layer_PEM(iref))/(inertiedat_PEM(ig,iref)**2))+ & ((layer_PEM(iref+1)-ice_table(ig))/(ice_inertia**2)))) else ! first layer is already a mixed layer ! (ie: take layer(iref=0)=0) therm_i(ig,1)=sqrt((layer_PEM(1))/ & (((ice_table(ig))/(inertiedat_PEM(ig,1)**2))+ & ((layer_PEM(1)-ice_table(ig))/(ice_inertia**2)))) endif ! of if (iref.ne.0) ! lower layers of pure ice do isoil=iref+2,nsoil_PEM therm_i(ig,isoil)=ice_inertia enddo endif ! of if (iref.lt.(nsoilmx)) endif ! permanent glaciers call soil_pem_1D(nsoil_PEM,.true.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) call soil_pem_1D(nsoil_PEM,.false.,therm_i(ig,:),timestep,tsurf_ave(ig),tsoil_ave(ig,:),alph_PEM,beta_PEM) do it = 1,timelen tsoil_inst(ig,:,it) = tsoil_inst(ig,:,it) - (Tcol_saved(:) - tsoil_ave(ig,:)) enddo error_depth = abs(icedepth_prev - ice_table(ig)) endif enddo error_depth = 1. countloop = 0 enddo END SUBROUTINE ini_icetable subroutine soil_pem_1D(nsoil,firstcall, & therm_i, & timestep,tsurf,tsoil,alph,beta) use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, & mu_PEM,fluxgeo use comsoil_h,only: volcapa implicit none !----------------------------------------------------------------------- ! Author: LL ! Purpose: Compute soil temperature using an implict 1st order scheme ! ! Note: depths of layers and mid-layers, soil thermal inertia and ! heat capacity are commons in comsoil_PEM.h ! A convergence loop is added until equilibrium !----------------------------------------------------------------------- #include "dimensions.h" !#include "dimphys.h" !#include"comsoil.h" !----------------------------------------------------------------------- ! arguments ! --------- ! inputs: integer,intent(in) :: nsoil ! number of soil layers logical,intent(in) :: firstcall ! identifier for initialization call real,intent(in) :: therm_i(nsoil) ! thermal inertia real,intent(in) :: timestep ! time step real,intent(in) :: tsurf ! surface temperature ! outputs: real,intent(inout) :: tsoil(nsoil) ! soil (mid-layer) temperature real,intent(inout) :: alph(nsoil-1) real,intent(inout) :: beta(nsoil-1) ! local variables: integer ik real :: thermdiff_PEM(nsoil-1) real :: mthermdiff_PEM(0:nsoil-1) real :: coefd_PEM(nsoil-1) real :: coefq_PEM(0:nsoil-1) ! 0. Initialisations and preprocessing step if (firstcall) then ! 0.1 Build mthermdiff_PEM(:), the mid-layer thermal diffusivities do ik=0,nsoil-1 mthermdiff_PEM(ik)=therm_i(ik+1)*therm_i(ik+1)/volcapa enddo ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities do ik=1,nsoil-1 thermdiff_PEM(ik)=((layer_PEM(ik)-mlayer_PEM(ik-1))*mthermdiff_PEM(ik) & +(mlayer_PEM(ik)-layer_PEM(ik))*mthermdiff_PEM(ik-1)) & /(mlayer_PEM(ik)-mlayer_PEM(ik-1)) enddo ! 0.3 Build coefficients mu_PEM, q_{k+1/2}, d_k, alpha_k and capcal ! mu_PEM mu_PEM=mlayer_PEM(0)/(mlayer_PEM(1)-mlayer_PEM(0)) ! q_{1/2} coefq_PEM(0)=volcapa*layer_PEM(1)/timestep ! q_{k+1/2} do ik=1,nsoil-1 coefq_PEM(ik)=volcapa*(layer_PEM(ik+1)-layer_PEM(ik)) & /timestep enddo ! d_k do ik=1,nsoil-1 coefd_PEM(ik)=thermdiff_PEM(ik)/(mlayer_PEM(ik)-mlayer_PEM(ik-1)) enddo ! alph_PEM_{N-1} alph(nsoil-1)=coefd_PEM(nsoil-1)/ & (coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1)) ! alph_PEM_k do ik=nsoil-2,1,-1 alph(ik)=coefd_PEM(ik)/(coefq_PEM(ik)+coefd_PEM(ik+1)* & (1.-alph(ik+1))+coefd_PEM(ik)) enddo endif ! of if (firstcall) IF (.not.firstcall) THEN ! 2. Compute soil temperatures ! First layer: tsoil(1)=(tsurf+mu_PEM*beta(1)* & thermdiff_PEM(1)/mthermdiff_PEM(0))/ & (1.+mu_PEM*(1.0-alph(1))*& thermdiff_PEM(1)/mthermdiff_PEM(0)) ! Other layers: do ik=1,nsoil-1 tsoil(ik+1)=alph(ik)*tsoil(ik)+beta(ik) enddo ENDIF ! 2. Compute beta_PEM coefficients (preprocessing for next time step) ! Bottom layer, beta_PEM_{N-1} beta(nsoil-1)=coefq_PEM(nsoil-1)*tsoil(nsoil) & /(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1)) & + fluxgeo/(coefq_PEM(nsoil-1)+coefd_PEM(nsoil-1)) ! Other layers do ik=nsoil-2,1,-1 beta(ik)=(coefq_PEM(ik)*tsoil(ik+1)+ & coefd_PEM(ik+1)*beta(ik+1))/ & (coefq_PEM(ik)+coefd_PEM(ik+1)*(1.0-alph(ik+1)) & +coefd_PEM(ik)) enddo end END MODULE ini_soil_mod