subroutine setgrid(nz,z,zmax,zfac) C construct regularly or geometrically spaced 1D grid C z(n)-z(1) = 2*z(1)*(zfac**(n-1)-1)/(zfac-1) implicit none integer NMAX parameter (NMAX=2000) integer nz, i real*8 zfac, zmax, z(NMAX), dz do i=1,nz dz = zmax/nz z(i) = (i-0.5)*dz enddo if (zfac>1.) then dz=zmax/(3.+2.*zfac*(zfac**(nz-2)-1.)/(zfac-1.)) z(1)=dz z(2)=3*z(1) do i=3,nz z(i)=(1.+zfac)*z(i-1)-zfac*z(i-2) enddo endif if (z(1)<1.e-5) then ! print *,'WARNING: first grid point is too shallow' endif end real*8 function smartzfac(nz_max,zmax,nz_delta,delta) C output can be used as input to setgrid C produces zfac with desired number of points within depth delta implicit none integer nz_max, nz_delta real*8 zmax, delta integer j real*8 f,g,gprime if (nz_max2.) then print *,'zfac=',smartzfac stop 'unwanted result in smartzfac' endif end subroutine smartgrid(nz,z,zdepth,thIn,rhoc,porosity,ti,rhocv, & layertype,icefrac) C*********************************************************************** C smartgrid: returns intelligently spaced grid and appropriate C values of thermal inertia ti and rho*c in icy layer C C INPUTS: C nz = number of grid points C z = grid spacing for dry regolith C (will be partially overwritten) C zdepth = depth (in meters) where ice table starts C negative values indicate no ice C rhoc = heat capacity per volume of dry regolith [J/m^3] C thIn = thermal inertia of dry regolith [SI-units] C porosity = void space / total volume C C OUTPUTS: z = grid coordinates C vectors ti and rhocv C C Earlier versions of this subroutine used layertype=1 C*********************************************************************** implicit none integer NMAX parameter (NMAX=2000) integer nz, j, b, layertype real*8 z(NMAX), zdepth, thIn, rhoc, porosity real*8 ti(NMAX), rhocv(NMAX), stretch, newrhoc, newti real*8 NULL, icefrac parameter (NULL=0.) if (zdepth>0.and.zdepth