Index: /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90 (revision 3770) @@ -19,24 +19,24 @@ ! Physical parameters -real, parameter :: d_dust = 5.78e-2 ! Tendency of dust [kg.m-2.y-1] -real, parameter :: dry_lag_porosity = 0.4 ! Porosity of dust lag +real, parameter :: d_dust = 5.78e-2 ! Tendency of dust [kg.m-2.y-1] +real, parameter :: dry_lag_porosity = 0.4 ! Porosity of dust lag real, parameter :: dry_regolith_porosity = 0.4 ! Porosity of regolith -real, parameter :: co2ice_porosity = 0 ! Porosity of CO2 ice -real, parameter :: h2oice_porosity = 0 ! Porosity of H2O ice -real, parameter :: rho_dust = 2500. ! Density of dust [kg.m-3] -real, parameter :: rho_rock = 3200. ! Density of rock [kg.m-3] +real, parameter :: breccia_porosity = 0.1 ! Porosity of breccia +real, parameter :: co2ice_porosity = 0. ! Porosity of CO2 ice +real, parameter :: h2oice_porosity = 0. ! Porosity of H2O ice +real, parameter :: rho_dust = 2500. ! Density of dust [kg.m-3] +real, parameter :: rho_rock = 3200. ! Density of rock [kg.m-3] ! Lag layer parameters -> see Levrard et al. 2007 -real, parameter :: hmin_lag = 0 ! Minimal height of the lag deposit to reduce the sublimation rate -real, parameter :: fred_subl = 1 ! Reduction factor of sublimation rate +real, parameter :: hmin_lag = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate +real, parameter :: fred_subl = 0.1 ! Reduction factor of sublimation rate ! Stratum = node of the linked list type :: stratum - real :: thickness ! Layer thickness [m] real :: top_elevation ! Layer top_elevation (top height from the surface) [m] - real :: co2ice_volfrac ! CO2 ice volumetric fraction - real :: h2oice_volfrac ! H2O ice volumetric fraction - real :: dust_volfrac ! Dust volumetric fraction - real :: pore_volfrac ! Pore volumetric fraction = pores + real :: h_co2ice ! Height of CO2 ice [m] + real :: h_h2oice ! Height of H2O ice [m] + real :: h_dust ! Height of dust [m] + real :: h_pore ! Height of pores [m] real :: poreice_volfrac ! Pore-filled ice volumetric fraction type(stratum), pointer :: up => null() ! Upper stratum (next node) @@ -70,10 +70,16 @@ ! > stratif2array ! > array2stratif -! > print_stratif +! > print_deposits +! Procedures to get information about a stratum: +! > thickness +! > is_co2ice_str +! > is_h2oice_str +! > is_dust_str +! > is_dust_lag +! > compute_h_toplag ! Procedures for the algorithm to evolve the layering: ! > make_layering -! > subl_co2ice_layering -! > subl_h2oice_layering -! > lift_dust_lags +! > lift_dust_lag +! > subl_ice_str !======================================================================= ! To display a stratum @@ -86,11 +92,10 @@ !---- Code -write(*,'(a,es13.6)') 'Thickness : ', this%thickness -write(*,'(a,es13.6)') 'Top elevation : ', this%top_elevation -write(*,'(a,es13.6)') 'CO2 ice (m3/m3) : ', this%co2ice_volfrac -write(*,'(a,es13.6)') 'H2O ice (m3/m3) : ', this%h2oice_volfrac -write(*,'(a,es13.6)') 'Dust (m3/m3) : ', this%dust_volfrac -write(*,'(a,es13.6)') 'Porosity (m3/m3) : ', this%pore_volfrac -write(*,'(a,es13.6)') 'Pore-filled ice (m3/m3): ', this%poreice_volfrac +write(*,'(a,es13.6)') 'Top elevation [m]: ', this%top_elevation +write(*,'(a,es13.6)') 'Height of CO2 ice [m]: ', this%h_co2ice +write(*,'(a,es13.6)') 'Height of H2O ice [m]: ', this%h_h2oice +write(*,'(a,es13.6)') 'Height of dust [m]: ', this%h_dust +write(*,'(a,es13.6)') 'Height of pores [m]: ', this%h_pore +write(*,'(a,es13.6)') 'Pore-filled ice [m3/m3]: ', this%poreice_volfrac END SUBROUTINE print_stratum @@ -98,5 +103,5 @@ !======================================================================= ! To add a stratum to the top of a layering -SUBROUTINE add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice) +SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice) implicit none @@ -104,5 +109,5 @@ !---- Arguments type(layering), intent(inout) :: this -real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice +real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice !---- Local variables @@ -113,24 +118,16 @@ allocate(str) nullify(str%up,str%down) -str%thickness = 1. str%top_elevation = 1. -str%co2ice_volfrac = 0. -str%h2oice_volfrac = 0. -str%dust_volfrac = 0. -str%pore_volfrac = 1. +str%h_co2ice = 0. +str%h_h2oice = 0. +str%h_dust = 1. +str%h_pore = 0. str%poreice_volfrac = 0. -if (present(thickness)) str%thickness = thickness if (present(top_elevation)) str%top_elevation = top_elevation -if (present(co2ice)) str%co2ice_volfrac = co2ice -if (present(h2oice)) str%h2oice_volfrac = h2oice -if (present(dust)) str%dust_volfrac = dust -if (present(pore)) str%pore_volfrac = pore +if (present(co2ice)) str%h_co2ice = co2ice +if (present(h2oice)) str%h_h2oice = h2oice +if (present(dust)) str%h_dust = dust +if (present(pore)) str%h_pore = pore if (present(poreice)) str%poreice_volfrac = poreice - -! Verification of volume fraction -if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then - call print_stratum(str) - error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!' -endif ! Increment the number of strata @@ -151,5 +148,5 @@ !======================================================================= ! To insert a stratum after another one in a layering -SUBROUTINE insert_stratum(this,str_prev,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice) +SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,poreice) implicit none @@ -158,5 +155,5 @@ type(layering), intent(inout) :: this type(stratum), pointer, intent(inout) :: str_prev -real, intent(in), optional :: thickness, top_elevation, co2ice, h2oice, dust, pore, poreice +real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice !---- Local variables @@ -165,29 +162,21 @@ !---- Code if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering - call add_stratum(this,thickness,top_elevation,co2ice,h2oice,dust,pore,poreice) + call add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice) else ! Creation of the stratum allocate(str) nullify(str%up,str%down) - str%thickness = 1. str%top_elevation = 1. - str%co2ice_volfrac = 0. - str%h2oice_volfrac = 0. - str%dust_volfrac = 0. - str%pore_volfrac = 1. + str%h_co2ice = 0. + str%h_h2oice = 0. + str%h_dust = 1. + str%h_pore = 0. str%poreice_volfrac = 0. - if (present(thickness)) str%thickness = thickness if (present(top_elevation)) str%top_elevation = top_elevation - if (present(co2ice)) str%co2ice_volfrac = co2ice - if (present(h2oice)) str%h2oice_volfrac = h2oice - if (present(dust)) str%dust_volfrac = dust - if (present(pore)) str%pore_volfrac = pore + if (present(co2ice)) str%h_co2ice = co2ice + if (present(h2oice)) str%h_h2oice = h2oice + if (present(dust)) str%h_dust = dust + if (present(pore)) str%h_pore = pore if (present(poreice)) str%poreice_volfrac = poreice - - ! Verification of volume fraction - if (abs(1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac + str%pore_volfrac)) > tol) then - call print_stratum(str) - error stop 'add_stratum: properties for the new stratum are not possible (sum of volumetric fraction /= 1)!' - endif ! Increment the number of strata @@ -203,5 +192,5 @@ current => str%up do while (associated(current)) - current%top_elevation = current%down%top_elevation + current%thickness + current%top_elevation = current%down%top_elevation + thickness(current) current => current%up enddo @@ -224,5 +213,5 @@ !---- Code -! Protect the 3 sub-surface strata from removing +! Protect the sub-surface strata from removing if (str%top_elevation <= 0.) return @@ -274,4 +263,6 @@ SUBROUTINE ini_layering(this) +use comsoil_h_PEM, only: soil_pem, nsoilmx_PEM, layer_PEM, index_breccia, index_bedrock + implicit none @@ -279,12 +270,31 @@ type(layering), intent(inout) :: this -!---- Code -! Creation of three strata at the bottom of the layering -! 1) Dry porous deep regolith -!call add_stratum(this,18.,-12.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity) -! 2) Ice-cemented regolith (ice table) -!call add_stratum(this,10.,-2.,0.,1. - dry_regolith_porosity,dry_regolith_porosity,0.) -! 3) Dry porous regolith -call add_stratum(this,2.,0.,0.,0.,dry_regolith_porosity,1. - dry_regolith_porosity) +!---- Local variables +integer :: i +real :: h_soil, h_pore, h_tot + +!---- Code +! Creation of strata at the bottom of the layering to describe the sub-surface +if (soil_pem) then + do i = nsoilmx_PEM,index_bedrock,-1 + h_soil = layer_PEM(i) - layer_PEM(i - 1) ! No porosity + call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,0.,0.) + enddo + do i = index_bedrock - 1,index_breccia,-1 + h_tot = layer_PEM(i) - layer_PEM(i - 1) + h_pore = h_tot*breccia_porosity + h_soil = h_tot - h_pore + call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.) + enddo + do i = index_breccia - 1,2,-1 + h_tot = layer_PEM(i) - layer_PEM(i - 1) + h_pore = h_tot*dry_regolith_porosity + h_soil = h_tot - h_pore + call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.) + enddo + h_pore = layer_PEM(1)*dry_regolith_porosity + h_soil = layer_PEM(1) - h_pore + call add_stratum(this,0.,0.,0.,h_soil,h_pore,0.) +endif END SUBROUTINE ini_layering @@ -387,11 +397,10 @@ k = 1 do while (associated(current)) - stratif_array(ig,islope,k,1) = current%thickness - stratif_array(ig,islope,k,2) = current%top_elevation - stratif_array(ig,islope,k,3) = current%co2ice_volfrac - stratif_array(ig,islope,k,4) = current%h2oice_volfrac - stratif_array(ig,islope,k,5) = current%dust_volfrac - stratif_array(ig,islope,k,6) = current%pore_volfrac - stratif_array(ig,islope,k,7) = current%poreice_volfrac + stratif_array(ig,islope,k,1) = current%top_elevation + stratif_array(ig,islope,k,2) = current%h_co2ice + stratif_array(ig,islope,k,3) = current%h_h2oice + stratif_array(ig,islope,k,4) = current%h_dust + stratif_array(ig,islope,k,5) = current%h_pore + stratif_array(ig,islope,k,6) = current%poreice_volfrac current => current%up k = k + 1 @@ -419,14 +428,6 @@ do islope = 1,nslope do ig = 1,ngrid - stratif(ig,islope)%bottom%thickness = stratif_array(ig,islope,1,1) - stratif(ig,islope)%bottom%top_elevation = stratif_array(ig,islope,1,2) - stratif(ig,islope)%bottom%co2ice_volfrac = stratif_array(ig,islope,1,3) - stratif(ig,islope)%bottom%h2oice_volfrac = stratif_array(ig,islope,1,4) - stratif(ig,islope)%bottom%dust_volfrac = stratif_array(ig,islope,1,5) - stratif(ig,islope)%bottom%pore_volfrac = stratif_array(ig,islope,1,6) - stratif(ig,islope)%bottom%poreice_volfrac = stratif_array(ig,islope,1,7) - do k = 2,size(stratif_array,3) - if (all(abs(stratif_array(ig,islope,k,:)) < eps)) exit - call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6),stratif_array(ig,islope,k,7)) + do k = 1,size(stratif_array,3) + call add_stratum(stratif(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6)) enddo enddo @@ -436,6 +437,6 @@ !======================================================================= -! To display a stratification (layerings) -SUBROUTINE print_stratif(this,ngrid,nslope) +! To display deposits +SUBROUTINE print_deposits(this,ngrid,nslope) implicit none @@ -459,9 +460,114 @@ enddo -END SUBROUTINE print_stratif - +END SUBROUTINE print_deposits + +!======================================================================= +!======================================================================= +! To compute the thickness of a stratum +FUNCTION thickness(this) RESULT(h_str) + +implicit none + +!---- Arguments +type(stratum), intent(in) :: this +real :: h_str + +!---- Code +h_str = this%h_co2ice + this%h_h2oice + this%h_dust + this%h_pore + +END FUNCTION thickness + +!======================================================================= +! To know whether the stratum can be considered as a CO2 ice layer +FUNCTION is_co2ice_str(str) RESULT(co2ice_str) + +implicit none + +!---- Arguments +type(stratum), intent(in) :: str +logical :: co2ice_str + +!---- Code +co2ice_str = .false. +if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true. + +END FUNCTION is_co2ice_str + +!======================================================================= +! To know whether the stratum can be considered as a H2O ice layer +FUNCTION is_h2oice_str(str) RESULT(h2oice_str) + +implicit none + +!---- Arguments +type(stratum), intent(in) :: str +logical :: h2oice_str + +!---- Code +h2oice_str = .false. +if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true. + +END FUNCTION is_h2oice_str + +!======================================================================= +! To know whether the stratum can be considered as a dust layer +FUNCTION is_dust_str(str) RESULT(dust_str) + +implicit none + +!---- Arguments +type(stratum), intent(in) :: str +logical :: dust_str + +!---- Code +dust_str = .false. +if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true. + +END FUNCTION is_dust_str + +!======================================================================= +! To know whether the stratum can be considered as a dust lag +FUNCTION is_dust_lag(str) RESULT(dust_lag) + +implicit none + +!---- Arguments +type(stratum), intent(in) :: str +logical :: dust_lag + +!---- Code +dust_lag = .false. +if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true. + +END FUNCTION is_dust_lag + +!======================================================================= +! To get the top lag height +FUNCTION thickness_toplag(this) RESULT(h_toplag) + +implicit none + +!---- Arguments +type(layering), intent(in) :: this +real :: h_toplag + +!---- Local variables +type(stratum), pointer :: current + +!---- Code +h_toplag = 0. +current => this%top +! Is the considered stratum a dust lag? +do while (current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol .and. associated(current)) + h_toplag = h_toplag + thickness(current) + current => current%down +enddo + +END FUNCTION thickness_toplag + +!======================================================================= !======================================================================= ! Layering algorithm -SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,zshift_surf,new_lag1,new_lag2,zlag,stopPEM,current1,current2) +SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current) implicit none @@ -469,11 +575,10 @@ !---- Arguments ! Inputs -real, intent(in) :: d_co2ice, d_h2oice, d_dust ! Outputs type(layering), intent(inout) :: this -type(stratum), pointer, intent(inout) :: current1, current2 -logical, intent(inout) :: new_str, new_lag1, new_lag2 -integer, intent(inout) :: stopPEM +type(stratum), pointer, intent(inout) :: current +real, intent(inout) :: d_co2ice, d_h2oice +logical, intent(inout) :: new_str, new_lag real, intent(out) :: zshift_surf, zlag @@ -481,6 +586,6 @@ real :: dh_co2ice, dh_h2oice, dh_dust real :: h_co2ice_old, h_h2oice_old, h_dust_old -real :: thickness, h2subl, h2subl_tot, h2lift, h2lift_tot -real :: h_toplag +real :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot +real :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag type(stratum), pointer :: currentloc @@ -489,235 +594,162 @@ dh_h2oice = d_h2oice/rho_h2oice dh_dust = d_dust/rho_dust -zshift_surf = 0. +zshift_surf = this%top%top_elevation zlag = 0. -if (dh_dust < 0.) then ! Dust lifting only - if (abs(dh_co2ice) > eps .or. abs(dh_h2oice) > eps) error stop 'Situation not managed: dust lifting + CO2/H2O ice condensation/sublimation!' - write(*,'(a)') ' Stratification -> Dust lifting' +! DUST SEDIMENTATION with possibly ice condensation +if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice .and. dh_co2ice >= 0. .and. dh_h2oice >= 0.) then + write(*,*) '> Layering: dust sedimentation' + h_pore = dh_dust*dry_regolith_porosity/(1. - dry_regolith_porosity) + h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore + ! New stratum at the top of the layering + if (new_str) then + call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.) + new_str = .false. + else + this%top%top_elevation = this%top%top_elevation + h_tot + this%top%h_co2ice = this%top%h_co2ice + dh_co2ice + this%top%h_h2oice = this%top%h_h2oice + dh_h2oice + this%top%h_dust = this%top%h_dust + dh_dust + this%top%h_pore = this%top%h_pore + h_pore + endif +!----------------------------------------------------------------------- +! DUST LIFTING +!(with possibly ice sublimation???) +!else if (dh_dust < dh_co2ice .and. dh_dust < dh_h2oice .and. dh_co2ice <= 0. .and. dh_h2oice <= 0.) then +else if (dh_dust < 0. .and. abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then + write(*,*) '> Layering: dust lifting' h2lift_tot = abs(dh_dust) - do while (h2lift_tot > 0. .and. associated(current1)) + do while (h2lift_tot > 0. .and. associated(current)) ! Is the considered stratum a dust lag? - if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we can lift dust - h_dust_old = current1%dust_volfrac*current1%thickness - ! How much dust can be lifted in the considered dust lag? - if (h2lift_tot - h_dust_old <= eps) then ! There is enough to lift - h2lift = h2lift_tot - h2lift_tot = 0. - call lift_dust_lags(this,current1,h2lift) - else ! Only a fraction can be lifted and so we move to the underlying stratum - h2lift = h_dust_old - h2lift_tot = h2lift_tot - h2lift - call lift_dust_lags(this,current1,h2lift) - current1 => current1%down - endif - else ! No, we stop - write(*,'(a)') ' Stratification -> Dust lifting: no dust lag!' - stopPEM = 8 + !if ((current%h_co2ice > eps .and. abs(dh_co2ice) < tol) .or. (current%h_h2oice > eps .and. abs(dh_h2oice) < tol)) then ! No, there is ice cementing the dust + if (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) then ! No, there is ice cementing the dust + dh_dust = 0. ! The tendency is set to 0 + exit + else ! Yes, we can lift dust + h2lift = min(h2lift_tot,current%h_dust) ! How much dust can be lifted? + h2lift_tot = h2lift_tot - h2lift + call lift_dust_lag(this,current,h2lift) + if (h2lift_tot > 0.) current => current%down ! Still dust to be lifted so we move to the underlying stratum + endif + enddo + if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0 +!----------------------------------------------------------------------- +! ICE CONDENSATION +else if ((dh_co2ice > dh_dust .or. dh_h2oice > dh_dust) .and. dh_dust >= 0.) then + write(*,*) '> Layering: ice condensation' + ! Which porosity is considered? + if (dh_co2ice > dh_h2oice) then ! CO2 ice is dominant + h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity) + else ! H2O ice is dominant + h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity) + endif + h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore + ! New stratum at the top of the layering + if (new_str) then + call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.) + new_str = .false. + else + this%top%top_elevation = this%top%top_elevation + h_tot + this%top%h_co2ice = this%top%h_co2ice + dh_co2ice + this%top%h_h2oice = this%top%h_h2oice + dh_h2oice + this%top%h_dust = this%top%h_dust + dh_dust + this%top%h_pore = this%top%h_pore + h_pore + endif +!----------------------------------------------------------------------- +! ICE SUBLIMATION +else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. abs(dh_dust) < tol) then + write(*,*) '> Layering: ice sublimation' + h2subl_co2ice_tot = abs(dh_co2ice) + h2subl_h2oice_tot = abs(dh_h2oice) + do while (h2subl_co2ice_tot > 0. .and. h2subl_h2oice_tot > 0. .and. associated(current)) + h_co2ice_old = current%h_co2ice + h_h2oice_old = current%h_h2oice + h_dust_old = current%h_dust + +!~ ! How much is CO2 ice sublimation inhibited by the top dust lag? +!~ h_toplag = 0. +!~ if (associated(current%up)) then ! If there is an upper stratum +!~ currentloc => current%up +!~ ! Is the considered stratum a dust lag? +!~ do while (.not. (current%h_co2ice > 0. .or. current%h_h2oice > 0. .or. current%poreice_volfrac > 0.) .and. associated(currentloc%up)) +!~ h_toplag = h_toplag + thickness(currentloc) +!~ currentloc => currentloc%up +!~ enddo +!~ endif +!~ h2subl_co2ice_tot = h2subl_co2ice_tot*fred_subl**(h_toplag/hmin_lag) + + ! Is there CO2 ice to sublimate? + h2subl_co2ice = 0. + if (h_co2ice_old > 0. .and. h2subl_co2ice_tot > 0.) then + h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old) + h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice + endif + ! Is there H2O ice to sublimate? + h2subl_h2oice = 0. + if (h_h2oice_old > 0. .and. h2subl_h2oice_tot > 0.) then + h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old) + h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice + endif + ! Sublimation + if (h2subl_co2ice > 0. .or. h2subl_h2oice >0.) call subl_ice_str(this,current,h2subl_co2ice,h2subl_h2oice,zlag,new_lag) + ! Still ice to be sublimated and no more ice in the current stratum so we move to the underlying stratum + if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol) then + current => current%down + new_lag = .true. + else exit endif enddo - if (h2lift_tot > 0.) then - write(*,'(a)') ' Stratification -> Not enough dust in the lag layers to complete the lifting!' - stopPEM = 8 - endif - zshift_surf = dh_dust + h2lift_tot - + if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0 + if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0 +!----------------------------------------------------------------------- +!~ ! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION +!~ else if (dh_co2ice < 0. .and. dh_h2oice > 0. .and. (abs(dh_dust) < abs(dh_co2ice) .or. abs(dh_dust) < abs(dh_h2oice))) then + +! TO DO???? + +!----------------------------------------------------------------------- +! NOTHING HAPPENS +else if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then + ! We do nothing +!----------------------------------------------------------------------- +! UNUSUAL (WEIRD) SITUATION else - -!------ Dust sedimentation only - if (abs(dh_co2ice) < eps .and. abs(dh_h2oice) < eps) then - write(*,'(a)') ' Stratification -> Dust sedimentation' - ! New stratum at the layering top by sedimentation of dust - thickness = dh_dust/(1. - dry_regolith_porosity) - if (thickness > 0.) then ! Only if the stratum is building up - if (new_str) then - call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,0.,dh_dust/thickness,dry_regolith_porosity,0.) - new_str = .false. - else - this%top%thickness = this%top%thickness + thickness - this%top%top_elevation = this%top%top_elevation + thickness - endif - endif - zshift_surf = thickness - -!------ Condensation of CO2 ice + H2O ice - else if ((dh_co2ice >= 0. .and. dh_h2oice > 0.) .or. (dh_co2ice > 0. .and. dh_h2oice >= 0.)) then - write(*,'(a)') ' Stratification -> CO2 and H2O ice condensation' - ! New stratum at the layering top by condensation of CO2 and H2O ice - thickness = dh_co2ice/(1. - co2ice_porosity) + dh_h2oice/(1. - h2oice_porosity) + dh_dust - if (thickness > 0.) then ! Only if the stratum is building up - if (new_str) then - call add_stratum(this,thickness,this%top%top_elevation + thickness,dh_co2ice/thickness,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_co2ice/thickness + dh_h2oice/thickness + dh_dust/thickness),0.) - new_str = .false. - else - this%top%thickness = this%top%thickness + thickness - this%top%top_elevation = this%top%top_elevation + thickness - endif - endif - zshift_surf = thickness - -!------ Sublimation of CO2 ice + Condensation of H2O ice - else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then - write(*,'(a)') ' Stratification -> Sublimation of CO2 ice + Condensation of H2O ice' - ! CO2 ice sublimation in the considered stratum + New stratum for dust lag - h2subl_tot = abs(dh_co2ice) - do while (h2subl_tot > 0. .and. associated(current1)) - h_co2ice_old = current1%co2ice_volfrac*current1%thickness - h_h2oice_old = current1%h2oice_volfrac*current1%thickness - h_dust_old = current1%dust_volfrac*current1%thickness - - ! How much is CO2 ice sublimation inhibited by the top dust lag? - h_toplag = 0. - if (associated(current1%up)) then ! If there is an upper stratum - currentloc => current1%up - do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol) - h_toplag = h_toplag + currentloc%thickness - if (.not. associated(currentloc%up)) exit - currentloc => currentloc%up - enddo - endif - h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag) - - ! How much CO2 ice can sublimate in the considered stratum? - if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate - h2subl = h2subl_tot - h2subl_tot = 0. - call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) - else if (h_co2ice_old < eps) then ! There is nothing to sublimate - ! Is the considered stratum a dust lag? - if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum - current1 => current1%down - new_lag1 = .true. - else ! No, so we stop here - exit - endif - else ! Only a fraction can sublimate and so we move to the underlying stratum - h2subl = h_co2ice_old - h2subl_tot = h2subl_tot - h2subl - call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) - current1 => current1%down - new_lag1 = .true. - endif - enddo - if (h2subl_tot > 0.) then - write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!' - stopPEM = 8 - endif - ! New stratum at the layering top by condensation of H2O ice - thickness = dh_h2oice/(1. - h2oice_porosity) + dh_dust - if (thickness > 0.) then ! Only if the stratum is building up - if (new_str) then - call add_stratum(this,thickness,this%top%top_elevation + thickness,0.,dh_h2oice/thickness,dh_dust/thickness,1. - (dh_h2oice/thickness + dh_dust/thickness),0.) - new_str = .false. - else - this%top%thickness = this%top%thickness + thickness - this%top%top_elevation = this%top%top_elevation + thickness - endif - endif - zshift_surf = zshift_surf + thickness - -!------ Sublimation of CO2 ice + H2O ice - else if ((dh_co2ice <= 0. .and. dh_h2oice < 0.) .or. (dh_co2ice < 0. .and. dh_h2oice <= 0.)) then - write(*,'(a)') ' Stratification -> Sublimation of CO2 and H2O ice' - ! CO2 ice sublimation in the considered stratum + New stratum for dust lag - h2subl_tot = abs(dh_co2ice) - do while (h2subl_tot > 0. .and. associated(current1)) - h_co2ice_old = current1%co2ice_volfrac*current1%thickness - h_h2oice_old = current1%h2oice_volfrac*current1%thickness - h_dust_old = current1%dust_volfrac*current1%thickness - - ! How much is CO2 ice sublimation inhibited by the top dust lag? - h_toplag = 0. - if (associated(current1%up)) then ! If there is an upper stratum - currentloc => current1%up - do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol) - h_toplag = h_toplag + currentloc%thickness - if (.not. associated(currentloc%up)) exit - currentloc => currentloc%up - enddo - endif - h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag) - - ! How much CO2 ice can sublimate in the considered stratum? - if (h2subl_tot - h_co2ice_old <= eps) then ! There is enough to sublimate - h2subl = h2subl_tot - h2subl_tot = 0. - call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) - else if (h_co2ice_old < eps) then ! There is nothing to sublimate - ! Is the considered stratum a dust lag? - if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum - current1 => current1%down - new_lag1 = .true. - else ! No, so we stop here - exit - endif - else ! Only a fraction can sublimate and so we move to the underlying stratum - h2subl = h_co2ice_old - h2subl_tot = h2subl_tot - h2subl - call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) - current1 => current1%down - new_lag1 = .true. - endif - enddo - if (h2subl_tot > 0.) then - write(*,'(a)') ' Stratification -> Not enough CO2 ice in the layering to complete the sublimation!' - stopPEM = 8 - endif - ! H2O ice sublimation in the considered stratum + New stratum for dust lag - h2subl_tot = abs(dh_h2oice) - do while (h2subl_tot > 0. .and. associated(current2)) - h_co2ice_old = current2%co2ice_volfrac*current2%thickness - h_h2oice_old = current2%h2oice_volfrac*current2%thickness - h_dust_old = current2%dust_volfrac*current2%thickness - - ! How much is H2O ice sublimation inhibited by the top dust lag? - h_toplag = 0. - if (associated(current2%up)) then ! If there is an upper stratum - currentloc => current2%up - do while (abs(1. - (currentloc%dust_volfrac + currentloc%pore_volfrac)) < tol) - h_toplag = h_toplag + currentloc%thickness - if (.not. associated(currentloc%up)) exit - currentloc => currentloc%up - enddo - endif - h2subl_tot = h2subl_tot*fred_subl**(h_toplag/hmin_lag) - - ! How much H2O ice can sublimate in the considered stratum? - if (h2subl_tot - h_h2oice_old <= eps) then ! There is enough to sublimate - h2subl = h2subl_tot - h2subl_tot = 0. - call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag) - else if (h_h2oice_old < eps) then ! There is nothing to sublimate - ! Is the considered stratum a dust lag? - if (abs(1. - (current1%dust_volfrac + current1%pore_volfrac)) < tol) then ! Yes, we move to the underlying stratum - current1 => current1%down - new_lag1 = .true. - else ! No, so we stop here - exit - endif - else ! Only a fraction can sublimate and so we move to the underlying stratum - h2subl = h_h2oice_old - h2subl_tot = h2subl_tot - h2subl - call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag) - current2 => current2%down - new_lag2 = .true. - endif - enddo - if (h2subl_tot > 0.) then - write(*,'(a)') ' Stratification -> Not enough H2O ice in the layering to complete the sublimation!' - stopPEM = 8 - endif - -!------ Condensation of CO2 ice + Sublimation of H2O ice - else if (dh_co2ice > 0. .and. dh_h2oice < 0.) then - error stop 'Impossible situation: CO2 ice condensation + H2O ice sublimation!' - endif + write(*,'(a)') 'Error: situation for the layering construction not managed!' + write(*,'(a,es12.3)') ' > dh_co2ice [m/y] = ', dh_co2ice + write(*,'(a,es12.3)') ' > dh_h2oice [m/y] = ', dh_h2oice + write(*,'(a,es12.3)') ' > dh_dust [m/y] = ', dh_dust, tol +!~ error stop endif +zshift_surf = this%top%top_elevation - zshift_surf + END SUBROUTINE make_layering !======================================================================= -! To sublimate CO2 ice in the layering -SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag) +! To lift dust in a stratum +SUBROUTINE lift_dust_lag(this,str,h2lift) + +implicit none + +!---- Arguments +type(layering), intent(inout) :: this +type(stratum), pointer, intent(inout) :: str +real, intent(in) :: h2lift + +!---- Code +! Update of properties in the eroding dust lag +str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust) +str%h_pore = str%h_pore - h2lift*str%h_pore/str%h_dust +str%h_dust = str%h_dust - h2lift + +! Remove the eroding dust lag if there is no dust anymore +if (str%h_dust < tol) call rm_stratum(this,str) + +END SUBROUTINE lift_dust_lag + +!======================================================================= +! To sublimate ice in a stratum +SUBROUTINE subl_ice_str(this,str,h2subl_co2ice,h2subl_h2oice,zlag,new_lag) implicit none @@ -727,131 +759,49 @@ type(stratum), pointer, intent(inout) :: str logical, intent(inout) :: new_lag -real, intent(inout) :: zshift_surf, zlag -real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old - -!---- Local variables -real :: r_subl, hsubl_dust, h_lag +real, intent(inout) :: zlag +real, intent(in) :: h2subl_co2ice, h2subl_h2oice + +!---- Local variables +real :: hsubl_dust, h_pore, h_ice, h_tot, h2subl type(stratum), pointer :: current !---- Code -! How much dust does CO2 ice sublimation involve to create the lag? -r_subl = h2subl/(1. - co2ice_porosity) & - /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity))) -hsubl_dust = r_subl*str%dust_volfrac*str%thickness - -str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust -zshift_surf = zshift_surf - h2subl/(1. - co2ice_porosity) - hsubl_dust -if (str%thickness < tol) then - ! Remove the sublimating stratum if there is no ice anymore - call rm_stratum(this,str) -else - ! Update of properties in the sublimating stratum - str%top_elevation = str%top_elevation - h2subl/(1. - co2ice_porosity) - hsubl_dust - str%co2ice_volfrac = (h_co2ice_old - h2subl)/str%thickness - str%h2oice_volfrac = h_h2oice_old/str%thickness - str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness - str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac) - str%poreice_volfrac = 0. - ! Correct the value of top_elevation for the upper strata - current => str%up - do while (associated(current)) - current%top_elevation = current%down%top_elevation + current%thickness - current => current%up - enddo -endif - -! New stratum = dust lag -h_lag = hsubl_dust/(1. - dry_lag_porosity) -if (h_lag > 0.) then ! Only if the dust lag is building up +h_ice = str%h_co2ice + str%h_h2oice +h2subl = h2subl_co2ice + h2subl_h2oice + +! How much dust does ice sublimation involve to create the lag? +hsubl_dust = str%h_dust*h2subl/h_ice + +! Update of properties in the sublimating stratum +str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hsubl_dust +str%h_co2ice = str%h_co2ice - h2subl_co2ice +str%h_h2oice = str%h_h2oice - h2subl_h2oice +str%h_dust = str%h_dust - hsubl_dust +str%h_pore = str%h_pore - h2subl*h_pore/h_ice + +! Correct the value of top_elevation for the upper strata +current => str%up +do while (associated(current)) + current%top_elevation = current%down%top_elevation + thickness(current) + current => current%up +enddo + +! Remove the sublimating stratum if there is no ice anymore +if (str%h_co2ice < tol .and. str%h_h2oice < tol) call rm_stratum(this,str) + +! New stratum above the current stratum as a dust lag +if (hsubl_dust > 0.) then ! Only if the dust lag is building up + h_pore = hsubl_dust*dry_lag_porosity/(1. - dry_lag_porosity) + h_tot = hsubl_dust + h_pore if (new_lag) then - call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity) + call insert_stratum(this,str,str%top_elevation + h_tot,0.,0.,hsubl_dust,h_pore) new_lag = .false. else - str%up%thickness = str%up%thickness + h_lag - str%up%top_elevation = str%up%top_elevation + h_lag + str%up%top_elevation = str%up%top_elevation + h_tot endif endif -zlag = zlag + h_lag - -END SUBROUTINE subl_co2ice_layering - -!======================================================================= -! To sublimate H2O ice in the layering -SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag) - -implicit none - -!---- Arguments -type(layering), intent(inout) :: this -type(stratum), pointer, intent(inout) :: str -logical, intent(inout) :: new_lag -real, intent(inout) :: zshift_surf, zlag -real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old - -!---- Local variables -real :: r_subl, hsubl_dust, h_lag -type(stratum), pointer :: current - -!---- Code -! How much dust does CO2 ice sublimation involve to create the lag? -r_subl = h2subl/(1. - h2oice_porosity) & - /(str%thickness*(str%co2ice_volfrac/(1. - co2ice_porosity) + str%h2oice_volfrac/(1. - h2oice_porosity))) -hsubl_dust = r_subl*str%dust_volfrac*str%thickness - -str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust -zshift_surf = zshift_surf - h2subl/(1. - h2oice_porosity) - hsubl_dust -if (str%thickness < tol) then - ! Remove the sublimating stratum if there is no ice anymore - call rm_stratum(this,str) -else - ! Update of properties in the sublimating stratum - str%top_elevation = str%top_elevation - h2subl/(1. - h2oice_porosity) - hsubl_dust - str%co2ice_volfrac = h_co2ice_old/str%thickness - str%h2oice_volfrac = (h_h2oice_old - h2subl)/str%thickness - str%dust_volfrac = (h_dust_old - hsubl_dust)/str%thickness - str%pore_volfrac = 1. - (str%co2ice_volfrac + str%h2oice_volfrac + str%dust_volfrac) - str%poreice_volfrac = 0. - ! Correct the value of top_elevation for the upper strata - current => str%up - do while (associated(current)) - current%top_elevation = current%down%top_elevation + current%thickness - current => current%up - enddo -endif - -! New stratum = dust lag -h_lag = hsubl_dust/(1. - dry_lag_porosity) -if (h_lag > 0.) then ! Only if the dust lag is building up - if (new_lag) then - call insert_stratum(this,str,h_lag,str%top_elevation + h_lag,0.,0.,1. - dry_lag_porosity,dry_lag_porosity) - new_lag = .false. - else - str%up%thickness = str%up%thickness + h_lag - str%up%top_elevation = str%up%top_elevation + h_lag - endif -endif -zlag = zlag + h_lag - -END SUBROUTINE subl_h2oice_layering - -!======================================================================= -! To lift dust in the layering -SUBROUTINE lift_dust_lags(this,str,h2lift) - -implicit none - -!---- Arguments -type(layering), intent(inout) :: this -type(stratum), pointer, intent(inout) :: str -real, intent(in) :: h2lift - -!---- Code -! Update of properties in the eroding dust lag -str%thickness = str%thickness - h2lift/(1. - str%pore_volfrac) -str%top_elevation = str%top_elevation - h2lift/(1. - str%pore_volfrac) -! Remove the eroding dust lag if there is no dust anymore -if (str%thickness < tol) call rm_stratum(this,str) - -END SUBROUTINE lift_dust_lags +zlag = zlag + h_tot + +END SUBROUTINE subl_ice_str END MODULE layering_mod Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90 (revision 3770) @@ -21,5 +21,5 @@ ! III Output -! III_a Update surface value for the PCM start files +! III_a Update surface values for the PCM start files ! III_b Write the "restart.nc" and "restartfi.nc" ! III_c Write the "restartpem.nc" @@ -69,5 +69,6 @@ use writediagpem_mod, only: writediagpem, writediagsoilpem use co2condens_mod, only: CO2cond_ps -use layering_mod, only: d_dust, ptrarray, stratum, layering, ini_layering, del_layering, make_layering, get_nb_str_max, nb_str_max, layering_algo +use layering_mod, only: ptrarray, stratum, layering, del_layering, make_layering, get_nb_str_max, & + nb_str_max, layering_algo, thickness_toplag, is_dust_lag, is_h2oice_str use dyn_ss_ice_m_mod, only: dyn_ss_ice_m use version_info_mod, only: print_version_info @@ -119,10 +120,7 @@ include "comgeom.h" include "iniprint.h" -include "callkeys.h" - -integer ngridmx -parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm) ! Same variable names as in the PCM +integer, parameter :: ngridmx = 2 + (jjm - 1)*iim - 1/jjm integer, parameter :: nlayer = llm ! Number of vertical layer integer :: ngrid ! Number of physical grid points @@ -196,8 +194,9 @@ real, dimension(:,:), allocatable :: q_co2_PEM_phys ! Grid points x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg] -! Variables for stratification (layering) evolution -type(layering), dimension(:,:), allocatable :: stratif ! Layering (linked list of "stratum") for each grid point and slope -type(ptrarray), dimension(:,:), allocatable :: current1, current2 ! Current active stratum in the layering -logical, dimension(:,:), allocatable :: new_str, new_lag1, new_lag2 ! Flags for the layering algorithm +! Variables for the evolution of layered deposits +type(layering), dimension(:,:), allocatable :: deposits ! Layering for each grid point and slope +type(ptrarray), dimension(:,:), allocatable :: current ! Current active stratum in the layering +logical, dimension(:,:), allocatable :: new_str, new_lag ! Flags for the layering algorithm +real :: h2o_ice_depth_old ! Old depth of subsurface ice layer ! Variables for slopes @@ -646,18 +645,10 @@ !------------------------ write(*,*) '> Reading "startpem.nc"' -allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),stratif(ngrid,nslope)) +allocate(co2_ice(ngrid,nslope),h2o_ice(ngrid,nslope),deposits(ngrid,nslope)) allocate(delta_h2o_adsorbed(ngrid),delta_co2_adsorbed(ngrid),delta_h2o_icetablesublim(ngrid)) -if (layering_algo) then - do islope = 1,nslope - do i = 1,ngrid - call ini_layering(stratif(i,islope)) - enddo - enddo -endif - delta_h2o_icetablesublim = 0. call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, & tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice, & - watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif) + watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,deposits) deallocate(tsurf_avg_yr1) @@ -667,8 +658,17 @@ co2ice_sublim_surf_ini = 0. h2oice_ini_surf = 0. +h2o_ice_depth = 0. is_co2ice_sublim_ini = .false. is_h2oice_sublim_ini = .false. is_co2ice_ini = .false. co2ice_disappeared = .false. +if (layering_algo) then + do ig = 1,ngrid + do islope = 1,nslope + co2_ice(ig,islope) = deposits(ig,islope)%top%h_co2ice + h2o_ice(ig,islope) = deposits(ig,islope)%top%h_h2oice + enddo + enddo +endif do i = 1,ngrid do islope = 1,nslope @@ -681,4 +681,5 @@ is_h2oice_sublim_ini(i,islope) = .true. h2oice_ini_surf = h2oice_ini_surf + cell_area(i)*subslope_dist(i,islope) + if (layering_algo) h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope)) endif enddo @@ -736,12 +737,10 @@ stopPEM = 0 if (layering_algo) then - allocate(new_str(ngrid,nslope),new_lag1(ngrid,nslope),new_lag2(ngrid,nslope),current1(ngrid,nslope),current2(ngrid,nslope)) + allocate(new_str(ngrid,nslope),new_lag(ngrid,nslope),current(ngrid,nslope)) new_str = .true. - new_lag1 = .true. - new_lag2 = .true. + new_lag = .true. do islope = 1,nslope do ig = 1,ngrid - current1(ig,islope)%p => stratif(ig,islope)%top - current2(ig,islope)%p => stratif(ig,islope)%top + current(ig,islope)%p => deposits(ig,islope)%top enddo enddo @@ -820,14 +819,16 @@ !------------------------ allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) - call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM) - call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) if (layering_algo) then do islope = 1,nslope do ig = 1,ngrid - call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),zshift_surf(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),zlag(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p) + call make_layering(deposits(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),new_str(ig,islope),zshift_surf(ig,islope),new_lag(ig,islope),zlag(ig,islope),current(ig,islope)%p) + co2_ice(ig,islope) = deposits(ig,islope)%top%h_co2ice + h2o_ice(ig,islope) = deposits(ig,islope)%top%h_h2oice enddo enddo else zlag = 0. + call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM) + call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) endif @@ -840,4 +841,12 @@ ps_timeseries,ps_avg_global_old,ps_avg_global_new,co2_ice,flag_co2flow) if (h2oice_flow .and. nslope > 1) call flow_h2oglaciers(ngrid,nslope,iflat,subslope_dist,def_slope_mean,tsurf_avg,h2o_ice,flag_h2oflow) + if (layering_algo) then + do islope = 1,nslope + do ig = 1,ngrid + deposits(ig,islope)%top%h_co2ice = co2_ice(ig,islope) + deposits(ig,islope)%top%h_h2oice = h2o_ice(ig,islope) + enddo + enddo + endif !------------------------ @@ -891,5 +900,5 @@ if (soil_pem) then ! II_d.2 Shifting soil temperature to surface - !call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM) + call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM) ! II_d.3 Update soil temperature @@ -1014,20 +1023,24 @@ call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,ps_avg_global_old,ps_avg_global_new) write(*,*) "> Updating the H2O sub-surface ice tendency due to lag layer" - do ig = 1,ngrid - do islope = 1,nslope - call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) - enddo - enddo + if (layering_algo) then + do ig = 1,ngrid + do islope = 1,nslope + if (is_h2oice_sublim_ini(ig,islope)) then + h2o_ice_depth_old = h2o_ice_depth(ig,islope) + h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope)) + call recomp_tend_h2o(h2o_ice_depth_old,h2o_ice_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) + endif + enddo + enddo +!~ else +!~ do ig = 1,ngrid +!~ do islope = 1,nslope +!~ call recomp_tend_h2o(icetable_depth_old(ig,islope),icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM_timeseries_old(ig,:,islope,:),tsoil_PEM_timeseries(ig,:,islope,:),d_h2oice(ig,islope)) +!~ enddo +!~ enddo + endif if (soil_pem) deallocate(icetable_depth_old,tsoil_PEM_timeseries_old) deallocate(vmr_co2_PEM_phys) write(*,*) "> Updating the H2O sub-surface ice depth" - do ig = 1,ngrid - do islope = 1,nslope - if (icetable_depth(ig,islope) > 0) then - icetable_depth(ig,islope)=icetable_depth(ig,islope) + d_h2oice(ig,islope)*0.01 !!! 0.01 is for 1 percent dust, needs to be updated !!! - endif - enddo - enddo - !------------------------ @@ -1036,6 +1049,4 @@ !------------------------ write(*,*) "> Checking the stopping criteria" - h2o_ice=sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) - co2_ice=sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid) call stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global_new,nslope) @@ -1062,6 +1073,4 @@ case(7) write(*,'(a,i0)') " **** STOPPING because the time limit for the PEM job will be reached soon: ", stopPEM - case(8) - write(*,'(a,i0)') " **** STOPPING because the layering algorithm met an hasty end: ", stopPEM case default write(*,'(a,i0)') " **** STOPPING with unknown stopping criterion code: ", stopPEM @@ -1080,5 +1089,5 @@ deallocate(d_co2ice,d_co2ice_ini,d_h2oice) deallocate(is_co2ice_ini,is_co2ice_sublim_ini,is_h2oice_sublim_ini) -if (layering_algo) deallocate(new_str,new_lag1,new_lag2,current1,current2) +if (layering_algo) deallocate(new_str,new_lag,current) !------------------------------ END RUN -------------------------------- @@ -1086,5 +1095,5 @@ !------------------------ ! III Output -! III_a Update surface value for the PCM start files +! III_a Update surface values for the PCM start files !------------------------ ! III_a.1 Ice update for start file @@ -1093,5 +1102,5 @@ write(*,*) '> Reconstructing perennial ice and frost for the PCM' watercap = 0. -perennial_co2ice =sum(stratif(ig,islope)%top%co2ice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) +perennial_co2ice = co2_ice do ig = 1,ngrid ! H2O ice metamorphism @@ -1102,15 +1111,11 @@ ! Is H2O ice still considered as an infinite reservoir for the PCM? - !if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then - if (sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then - ! There is enough ice to be considered as an infinite reservoir + if (sum(h2o_ice(ig,:)*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) > inf_h2oice_threshold) then + ! There is enough ice to be considered as an infinite reservoir watercaptag(ig) = .true. else - ! There too little ice to be considered as an infinite reservoir so ice is transferred to the frost + ! Too little ice to be considered as an infinite reservoir so ice is transferred to the frost watercaptag(ig) = .false. - - ! qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:) - - qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + sum(stratif(ig,islope)%top%h2oice_volfrac*stratif(ig,islope)%top%thickness*subslope_dist(ig,:)/cos(pi*def_slope_mean(:)/180.)) + qsurf(ig,igcm_h2o_ice,:) = qsurf(ig,igcm_h2o_ice,:) + h2o_ice(ig,:) h2o_ice(ig,:) = 0. endif @@ -1123,10 +1128,25 @@ enddo -! III.a.2. Tsurf update for start file +! III_a.2 Subsurface-surface interaction update for start file +zdqsdif_ssi_tot = 0. +h2o_ice_depth = 0. +if (layering_algo) then + write(*,*) '> Reconstructing ice sublimation flux from subsurface to surface for the PCM' + do ig = 1,ngrid + do islope = 1,nslope + if (is_dust_lag(deposits(ig,islope)%top) .and. is_h2oice_str(deposits(ig,islope)%top%down)) then + zdqsdif_ssi_tot(ig,islope) = d_h2oice(ig,islope) + h2o_ice_depth(ig,islope) = thickness_toplag(deposits(ig,islope)) + endif + enddo + enddo +endif + +! III.a.3. Tsurf update for start file write(*,*) '> Reconstructing the surface temperature for the PCM' tsurf = tsurf_avg + tsurf_dev deallocate(tsurf_dev) -! III_a.3 Tsoil update for start file +! III_a.4 Tsoil update for start file if (soil_pem) then write(*,*) '> Reconstructing the soil temperature profile for the PCM' @@ -1143,5 +1163,5 @@ deallocate(tsurf_avg,tsoil_dev) -! III_a.4 Pressure update for start file +! III_a.5 Pressure update for start file write(*,*) '> Reconstructing the pressure for the PCM' allocate(ps_start(ngrid)) @@ -1150,5 +1170,5 @@ deallocate(ps_avg,ps_dev) -! III_a.5 Tracers update for start file +! III_a.6 Tracers update for start file write(*,*) '> Reconstructing the tracer VMR for the PCM' allocate(zplev_start0(ngrid,nlayer + 1),zplev_new(ngrid,nlayer + 1)) @@ -1187,5 +1207,5 @@ q(ig,l + 1,nnq) = q(ig,l + 1,nnq) + extra_mass*(zplev_new(ig,l + 1) - zplev_new(ig,l + 2)) write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq) /= "dust_number",noms(nnq) /= "ccn_number" - endif + endif if (q(ig,l,nnq) < 0) q(ig,l,nnq) = 1.e-30 enddo @@ -1194,5 +1214,5 @@ deallocate(zplev_new) -! III_a.6 Albedo update for start file +! III_a.7 Albedo update for start file write(*,*) '> Reconstructing the albedo for the PCM' do ig = 1,ngrid @@ -1231,5 +1251,5 @@ enddo -! III_a.7 Orbital parameters update for start file +! III_a.8 Orbital parameters update for start file write(*,*) '> Setting the new orbital parameters' if (evol_orbit_pem) call recomp_orb_param(i_myear,i_myear_leg) @@ -1293,8 +1313,8 @@ !------------------------ write(*,*) '> Writing "restartpem.nc"' -if (layering_algo) nb_str_max = get_nb_str_max(stratif,ngrid,nslope) ! Get the maximum number of "stratum" in the stratification (layerings) +if (layering_algo) nb_str_max = get_nb_str_max(deposits,ngrid,nslope) ! Get the maximum number of "stratum" in the depositsication (layerings) call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & - co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif) + co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,deposits) call info_PEM(i_myear_leg,stopPEM,i_myear,n_myear) @@ -1302,7 +1322,7 @@ write(*,*) write(*,*) '****** PEM final information *******' -write(*,'(a,f10.2,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years." -write(*,'(a,f10.2,a,f10.2,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years." -write(*,'(a,f10.2,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years." +write(*,'(a,f16.4,a)') " + The PEM leg has run for ", i_myear_leg, " Martian years." +write(*,'(a,f16.4,a,f16.4,a)') " + The chained simulation has run for ", i_myear, " Martian years =", i_myear*convert_years, " Earth years." +write(*,'(a,f16.4,a)') " + The reached date is now ", (year_bp_ini + i_myear)*convert_years, " Earth years." write(*,*) "+ PEM: so far, so good!" write(*,*) '************************************' @@ -1311,10 +1331,10 @@ do islope = 1,nslope do i = 1,ngrid - call del_layering(stratif(i,islope)) + call del_layering(deposits(i,islope)) enddo enddo endif deallocate(q,longitude,latitude,cell_area,tsoil_PEM) -deallocate(co2_ice,h2o_ice,stratif) +deallocate(co2_ice,h2o_ice,deposits) !----------------------------- END OUTPUT ------------------------------ Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90 (revision 3770) @@ -21,5 +21,6 @@ use abort_pem_mod, only: abort_pem use compute_soiltemp_mod, only: ini_tsoil_pem, compute_tsoil_pem -use layering_mod, only: layering, array2stratif, nb_str_max, layering_algo +use layering_mod, only: layering, array2stratif, nb_str_max, layering_algo, add_stratum, ini_layering +use callkeys_mod, only: startphy_file #ifndef CPP_STD @@ -31,6 +32,4 @@ implicit none - -include "callkeys.h" character(*), intent(in) :: filename ! name of the startfi_PEM.nc @@ -134,35 +133,67 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Stratification (layerings) - nb_str_max = 1 + nb_str_max = 68 if (layering_algo) then found = inquire_dimension("nb_str_max") if (.not. found) then write(*,*) 'Pemetat0: failed loading !' - write(*,*) '''nb_str_max'' is set to 1!' + write(*,*) '''nb_str_max'' is set to 68 (sub-surface layers)!' else nb_str_max = int(inquire_dimension_length('nb_str_max')) endif - allocate(stratif_array(ngrid,nslope,nb_str_max,7)) + allocate(stratif_array(ngrid,nslope,nb_str_max,6)) stratif_array = 0. do islope = 1,nslope write(num,'(i2.2)') islope - call get_field('stratif_slope'//num//'_thickness',stratif_array(:,islope,:,1),found) - call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,2),found2) - found = found .or. found2 - call get_field('stratif_slope'//num//'_co2ice_volfrac',stratif_array(:,islope,:,3),found2) - found = found .or. found2 - call get_field('stratif_slope'//num//'_h2oice_volfrac',stratif_array(:,islope,:,4),found2) - found = found .or. found2 - call get_field('stratif_slope'//num//'_dust_volfrac',stratif_array(:,islope,:,5),found2) - found = found .or. found2 - call get_field('stratif_slope'//num//'_pore_volfrac',stratif_array(:,islope,:,6),found2) - found = found .or. found2 - call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,7),found2) - found = found .or. found2 - if (.not. found) then - write(*,*) 'Pemetat0: failed loading at least one of the properties of ' - endif ! not found + call get_field('stratif_slope'//num//'_top_elevation',stratif_array(:,islope,:,1),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' + exit + endif + call get_field('stratif_slope'//num//'_h_co2ice',stratif_array(:,islope,:,2),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' + exit + endif + call get_field('stratif_slope'//num//'_h_h2oice',stratif_array(:,islope,:,3),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' + exit + endif + call get_field('stratif_slope'//num//'_h_dust',stratif_array(:,islope,:,4),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' + exit + endif + call get_field('stratif_slope'//num//'_h_pore',stratif_array(:,islope,:,5),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' + exit + endif + call get_field('stratif_slope'//num//'_icepore_volfrac',stratif_array(:,islope,:,6),found) + if (.not. found) then + write(*,*) 'Pemetat0: failed loading !' + exit + endif enddo ! islope - call array2stratif(stratif_array,ngrid,nslope,stratif) + if (.not. found) then + write(*,*) 'So the deposits are initialized with sub-surface strata.' + write(*,*) 'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.' + do ig = 1,ngrid + if (watercaptag(ig)) then + do islope = 1,nslope + call ini_layering(stratif(ig,islope)) + call add_stratum(stratif(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.) + enddo + else + do islope = 1,nslope + call ini_layering(stratif(ig,islope)) + if (perennial_co2ice(ig,islope) > 0.) call add_stratum(stratif(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.) + enddo + endif + enddo + else + call array2stratif(stratif_array,ngrid,nslope,stratif) + endif deallocate(stratif_array) endif @@ -337,5 +368,5 @@ call close_startphy -else !No startpem, let's build all by hand +else ! No startpem, let's build all by hand write(*,*)'There is no "startpem.nc"!' @@ -349,10 +380,26 @@ ! co2 ice - write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice''.' + write(*,*)'So ''co2_ice'' is initialized with ''perennial_co2ice'' found in the PCM.' co2_ice = perennial_co2ice ! Stratification (layerings) - nb_str_max = 1 - if (layering_algo) write(*,*)'So the stratification (layerings) is initialized with 1 sub-surface stratum.' + nb_str_max = 68 + if (layering_algo) then + write(*,*)'So the deposits are initialized with sub-surface strata.' + write(*,*)'Ice is added with ''ini_huge_h2oice'' where ''watercaptag'' is true and otherwise with ''perennial_co2ice'' found in the PCM.' + do ig = 1,ngrid + if (watercaptag(ig)) then + do islope = 1,nslope + call ini_layering(stratif(ig,islope)) + call add_stratum(stratif(ig,islope),ini_huge_h2oice,0.,ini_huge_h2oice,0.,0.,0.) + enddo + else + do islope = 1,nslope + call ini_layering(stratif(ig,islope)) + if (perennial_co2ice(ig,islope) > 0.) call add_stratum(stratif(ig,islope),perennial_co2ice(ig,islope),perennial_co2ice(ig,islope),0.,0.,0.,0.) + enddo + endif + enddo + endif if (soil_pem) then Index: /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/pemredem.F90 (revision 3770) @@ -68,8 +68,4 @@ implicit none -#ifndef CPP_STD - include "callkeys.h" -#endif - character(*), intent(in) :: filename integer, intent(in) :: nsoil_PEM, ngrid, nslope @@ -99,15 +95,14 @@ if (layering_algo) then - allocate(stratif_array(ngrid,nslope,nb_str_max,7)) + allocate(stratif_array(ngrid,nslope,nb_str_max,6)) call stratif2array(stratif,ngrid,nslope,stratif_array) do islope = 1,nslope write(num,fmt='(i2.2)') islope - call put_field('stratif_slope'//num//'_thickness','Layering thickness',stratif_array(:,islope,:,1),Year) - call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,2),Year) - call put_field('stratif_slope'//num//'_co2ice_volfrac','Layering CO2 ice volume fraction',stratif_array(:,islope,:,3),Year) - call put_field('stratif_slope'//num//'_h2oice_volfrac','Layering H2O ice volume fraction',stratif_array(:,islope,:,4),Year) - call put_field('stratif_slope'//num//'_dust_volfrac','Layering dust volume fraction',stratif_array(:,islope,:,5),Year) - call put_field('stratif_slope'//num//'_pore_volfrac','Layering pore volume fraction',stratif_array(:,islope,:,6),Year) - call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,7),Year) + call put_field('stratif_slope'//num//'_top_elevation','Layering top elevation',stratif_array(:,islope,:,1),Year) + call put_field('stratif_slope'//num//'_h_co2ice','Layering CO2 ice height',stratif_array(:,islope,:,2),Year) + call put_field('stratif_slope'//num//'_h_h2oice','Layering H2O ice height',stratif_array(:,islope,:,3),Year) + call put_field('stratif_slope'//num//'_h_dust','Layering dust height',stratif_array(:,islope,:,4),Year) + call put_field('stratif_slope'//num//'_h_pore','Layering pore height',stratif_array(:,islope,:,5),Year) + call put_field('stratif_slope'//num//'_icepore_volfrac','Layering ice pore volume fraction',stratif_array(:,islope,:,6),Year) enddo deallocate(stratif_array) Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_mod.F90 (revision 3770) @@ -92,5 +92,9 @@ Rz_old = icetable_depth_old*zcdv/coef_diff ! Old resistance from PCM Rz_new = icetable_depth_new*zcdv/coef_diff ! New resistance based on new depth -R_dec = Rz_new/Rz_old ! Decrease because of resistance +if (abs(Rz_old) < 1.e-10) then + R_dec = 1. +else + R_dec = Rz_new/Rz_old ! Decrease because of resistance +endif ! The maxmimum of the daily averages over one year for the saturation vapor pressure at the ice table location @@ -103,7 +107,9 @@ ! Lower humidity due to growing lag layer (higher depth) +print*, psv_max_old,psv_max_new hum_dec = psv_max_old/psv_max_new ! Decrease because of lower water vapor pressure at the new depth ! Flux correction (decrease) +print*, 'aaaaah', R_dec,hum_dec d_h2oice = d_h2oice/R_dec/hum_dec Index: /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/soil_settings_PEM_mod.F90 (revision 3770) @@ -57,5 +57,5 @@ do iloop = 0,nsoil_PEM - 1 - if(lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then + if (lay1*(alpha**(iloop-0.5)) > mlayer_PEM(nsoil_PCM - 1)) then index_powerlaw = iloop exit @@ -69,5 +69,5 @@ ! Build layer_PEM() do iloop = 1,nsoil_PEM - 1 -layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2. + layer_PEM(iloop) = (mlayer_PEM(iloop) + mlayer_PEM(iloop - 1))/2. enddo layer_PEM(nsoil_PEM) = 2*mlayer_PEM(nsoil_PEM - 1) - mlayer_PEM(nsoil_PEM - 2) Index: /trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90 =================================================================== --- /trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90 (revision 3769) +++ /trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90 (revision 3770) @@ -69,6 +69,4 @@ endif -!if (abs(h2oice_ini_surf) < 1.e-5) stopPEM = 0 - END SUBROUTINE stopping_crit_h2o_ice @@ -133,6 +131,4 @@ endif -!if (abs(co2ice_sublim_surf_ini) < 1.e-5) stopPEM = 0 - if (ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)) then stopPEM = 4