!***************************************************** ! Commonly used subroutines for fast method ! written by Norbert Schorghofer 2007-2011 !***************************************************** pure function zint(y1,y2,z1,z2) ! interpolate between two values, y1*y2<0 implicit none real(8), intent(IN) :: y1,y2,z1,z2 real(8) zint zint = (y1*z2 - y2*z1)/(y1-y2) end function zint pure function colint(y,z,nz,i1,i2) ! column integrates y implicit none integer, intent(IN) :: nz, i1, i2 real(8), intent(IN) :: y(nz),z(nz) real(8) colint integer i real(8) dz(nz) dz(1)=z(1) ! this leaves z(1)/2 for surface point do i=2,nz-1 dz(i) = (z(i+1)-z(i-1))/2. enddo dz(nz) = z(nz)-z(nz-1) colint= sum(y(i1:i2)*dz(i1:i2)) end function colint pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1) !*********************************************************************** ! returns equilibrium depth for given ice content ! this is not the true (final) equilibrium depth !*********************************************************************** use allinterfaces, except_this_one => equildepth implicit none integer, intent(IN) :: nz real(8), intent(IN) :: z(nz), rhosatav(nz) real(8), intent(IN) :: rhosatav0, avrho1 integer i, typeE real(8) equildepth ! =zdepthE !real(8), external :: zint ! defined in allinterfaces.mod typeE = -9; equildepth = -9999. do i=1,nz if (rhosatav(i) <= avrho1) then typeE=i exit endif enddo if (typeE>1) then ! interpolate equildepth=zint(avrho1-rhosatav(typeE-1),avrho1-rhosatav(typeE),z(typeE-1),z(typeE)) end if if (typeE==1) then equildepth=zint(avrho1-rhosatav0,avrho1-rhosatav(1),0.d0,z(1)) if (equildepth>z(nz)) equildepth=-9999. endif end function equildepth subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1, & & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG) !*********************************************************************** ! returns interface depth and ypp ! also returns lower geothermally limited boundary, if applicable ! ! this is a crucial subroutine ! ! B = Diff*bigstep/(porosity*icedensity) [SI units] !*********************************************************************** use allinterfaces, except_this_one => depths_avmeth implicit none integer, intent(IN) :: nz, typeT real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill real(8), intent(IN) :: rhosatav0, rlow, avrho1 integer, intent(OUT) :: typeF ! index of depth below which filling occurs real(8), intent(INOUT) :: zdepthF real(8), intent(IN) :: B real(8), intent(OUT) :: ypp(nz), zdepthG integer, intent(INOUT) :: typeG ! positive on input when Fgeotherm>0 integer i, typeP, nlb, newtypeG real(8) eta(nz), Jpump1, help(nz), yp(nz), zdepthFold, ap_one, ap(nz) real(8) leak, cumfill, cumfillabove if (typeT<0) then nlb = nz do i=1,nz eta(i) = constriction(porefill(i)) enddo else !nlb = typeT-1 nlb = typeT ! changed 2010-09-29 do i=1,typeT-1 eta(i) = constriction(porefill(i)) enddo do i=typeT,nz eta(i)=0. enddo end if !-fill depth zdepthFold = zdepthF typeF = -9; zdepthF = -9999. call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp) ! yp also used below do i=1,nlb Jpump1 = (rhosatav(i)-avrho1)/z(i) ! <0 when stable ! yp is always <0 help(i) = Jpump1 - eta(i)*yp(i) leak = porefill(i)/B*(z(i)-zdepthFold)/(18./8314.) !help(i) = help(i)-leak ! optional if (help(i) <= 0.) then typeF=i !print *,'#',typeF,Jpump1,eta(typeF)*yp(typeF),leak exit endif enddo if (typeF>1) zdepthF = zint(help(typeF-1),help(typeF),z(typeF-1),z(typeF)) if (typeF==1) zdepthF=z(1) !-depth to shallowest perennial ice typeP = -9 do i=1,nz if (porefill(i)>0.) then typeP = i ! first point with ice exit endif enddo !-calculate ypp !call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp) call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap) if (typeP>0 .and. typeP0) typeG=-9 if (typeG<0) zdepthG=-9999. if (typeG>0 .and. typeT<0) then typeG=-9 do i=2,nz if (yp(i)>0.) then ! first point with reversed flux typeG=i zdepthG=zint(yp(i-1),yp(i),z(i-1),z(i)) !zdepthG=z(typeG) exit endif enddo else typeG = -9 endif if (typeG>0 .and. typeT<0) then cumfillabove = colint(porefill(:)/eta(:),z,nz,typeG-1,nz) newtypeG = -9 do i=typeG,nz if (minval(eta(i:nz))<=0.) cycle ! eta=0 means completely full cumfill=colint(porefill(:)/eta(:),z,nz,i,nz) if (cumfilltypeG) then write(34,*) '# adjustment to geotherm depth by',i-typeG zdepthG = zint(yp(i-1)*18./8314.*B-cumfillabove, & ! no good & yp(i)*18./8314.*B-cumfill,z(i-1),z(i)) if (zdepthG>z(i) .or. zdepthG0) typeG=newtypeG end if ! if typeG>0, then all ice at and below typeG should be erased end subroutine depths_avmeth pure function constriction(porefill) ! specify constriction function here, 0<=eta<=1 implicit none real(8), intent(IN) :: porefill real(8) eta, constriction if (porefill<=0.) eta = 1. if (porefill>0. .and. porefill<1.) then ! eta = 1. ! eta = 1-porefill eta = (1-porefill)**2 ! Hudson et al., JGR, 2009 endif if (porefill>=1.) eta = 0. constriction = eta end function constriction pure subroutine icechanges_poreonly(nz,z,typeF,typeG,avdrhoP,ypp,B,porefill) use allinterfaces implicit none integer, intent(IN) :: nz, typeF, typeG real(8), intent(IN) :: z(nz), ypp(nz), avdrhoP, B real(8), intent(INOUT) :: porefill(nz) integer j, erase, newtypeP, ub real(8) integ !----retreat ! avdrhoP>0 is outward loss from zdepthP ! avdrhoP<0 means gain at zdepthP or no ice anywhere if (avdrhoP>0.) then erase=0 do j=1,nz if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j) erase = j if (integ>B*avdrhoP*18./8314.) exit end do if (erase>0) porefill(1:erase)=0. endif ! new depth newtypeP = -9 do j=1,nz if (porefill(j)>0.) then newtypeP = j ! first point with ice exit endif enddo !----diffusive filling ub = typeF if (newtypeP>0 .and. typeF>0 .and. newtypeP0) then do j=ub,nz ! B=Diff/(porosity*icedensity)*86400*365.24*bigstep porefill(j) = porefill(j) + B*ypp(j) if (porefill(j)<0.) porefill(j)=0. if (porefill(j)>1.) porefill(j)=1. enddo end if !----enact bottom boundary if (typeG>0) porefill(typeG:nz)=0. end subroutine icechanges_poreonly pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, & & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG) !*********************************************************** ! advances ice table, advances interface, grows pore ice ! ! a crucial subroutine !*********************************************************** use miscparameters, only : icedensity use allinterfaces implicit none integer, intent(IN) :: nz, typeF, typeG real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP real(8), intent(IN) :: Diff, porosity, icefrac, bigstep real(8), intent(INOUT) :: zdepthT, porefill(nz) integer j, erase, newtypeP, ub, typeP, typeT real(8) B, beta, integ B = Diff*bigstep*86400.*365.24/(porosity*icedensity) ! advance ice table, avdrho>0 is retreat if (zdepthT>=0. .and. avdrho>0.) then typeP=-9999; typeT=-9999 do j=1,nz if (z(j)>zdepthT) then typeT=j exit endif enddo do j=1,nz if (porefill(j)>0.) then typeP=j exit endif enddo if (typeP==typeT) then ! new 2011-09-01 beta = 2*Diff*bigstep*(1-icefrac)/(1-porosity)*86400*365.24/icedensity zdepthT = sqrt(beta*avdrho*18./8314. + zdepthT**2) endif endif if (zdepthT>z(nz)) zdepthT=-9999. ! advance interface, avdrhoP>0 is loss from zdepthP if (avdrhoP>0.) then erase=0 do j=1,nz if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF if (zdepthT>=0. .and. z(j)>zdepthT) exit integ = colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j) erase = j if (integ>B*avdrhoP*18./8314.) exit end do if (erase>0) porefill(1:erase)=0. endif ! new depth newtypeP = -9 do j=1,nz if (zdepthT>=0. .and. z(j)>zdepthT) exit if (porefill(j)>0.) then newtypeP = j ! first point with pore ice exit endif enddo ! diffusive filling ub = typeF if (newtypeP>0 .and. typeF>0 .and. newtypeP0) then do j=ub,nz porefill(j) = porefill(j) + B*ypp(j) if (porefill(j)<0.) porefill(j)=0. if (porefill(j)>1.) porefill(j)=1. if (zdepthT>=0. .and. z(j)>zdepthT) exit enddo end if ! below icetable if (zdepthT>=0.) then do j=1,nz if (z(j)>zdepthT) porefill(j) = icefrac/porosity enddo else ! geothermal lower boundary if (typeG>0) porefill(typeG:nz)=0. end if end subroutine icechanges subroutine assignthermalproperties(nz,thIn,rhoc, & & ti,rhocv,typeT,icefrac,porosity,porefill) !********************************************************* ! assign thermal properties of soil !********************************************************* implicit none integer, intent(IN) :: nz integer, intent(IN), optional :: typeT real(8), intent(IN), optional :: icefrac real(8), intent(IN) :: thIn, rhoc real(8), intent(IN), optional :: porosity, porefill(nz) real(8), intent(OUT) :: ti(nz), rhocv(nz) integer j real(8) newrhoc, newti, fill real(8), parameter :: NULL=0. ti(1:nz) = thIn rhocv(1:nz) = rhoc if (typeT>0) then call soilthprop(porosity,NULL,rhoc,thIn,2,newrhoc,newti,icefrac) rhocv(typeT:nz) = newrhoc ti(typeT:nz) = newti endif do j=1,nz fill = porefill(j) ! off by half point if (fill>0. .and. (typeT<0 .or. (typeT>0 .and. j