| 1 | subroutine setgrid(nz,z,zmax,zfac) |
|---|
| 2 | ! construct regularly or geometrically spaced 1D grid |
|---|
| 3 | ! z(n)-z(1) = 2*z(1)*(zfac**(n-1)-1)/(zfac-1) |
|---|
| 4 | ! choice of z(1) and z(2) is compatible with conductionQ |
|---|
| 5 | implicit none |
|---|
| 6 | integer, intent(IN) :: nz |
|---|
| 7 | real(8), intent(IN) :: zfac, zmax |
|---|
| 8 | real(8), intent(OUT) :: z(nz) |
|---|
| 9 | integer i |
|---|
| 10 | real(8) dz |
|---|
| 11 | |
|---|
| 12 | dz = zmax/nz |
|---|
| 13 | do i=1,nz |
|---|
| 14 | z(i) = (i-0.5)*dz |
|---|
| 15 | enddo |
|---|
| 16 | |
|---|
| 17 | if (zfac>1.) then |
|---|
| 18 | dz = zmax/(3.+2.*zfac*(zfac**(nz-2)-1.)/(zfac-1.)) |
|---|
| 19 | z(1) = dz |
|---|
| 20 | z(2) = 3*z(1) |
|---|
| 21 | do i=3,nz |
|---|
| 22 | z(i) = (1.+zfac)*z(i-1) - zfac*z(i-2) |
|---|
| 23 | ! z(i) = z(i-1) + zfac*(z(i-1)-z(i-2)) ! equivalent |
|---|
| 24 | enddo |
|---|
| 25 | endif |
|---|
| 26 | |
|---|
| 27 | end subroutine setgrid |
|---|
| 28 | |
|---|
| 29 | |
|---|
| 30 | |
|---|
| 31 | real(8) function smartzfac(nz_max,zmax,nz_delta,delta) |
|---|
| 32 | ! output can be used as input to setgrid |
|---|
| 33 | ! produces zfac with desired number of points within depth delta |
|---|
| 34 | implicit none |
|---|
| 35 | integer, intent(IN) :: nz_max, nz_delta |
|---|
| 36 | real(8), intent(IN) :: zmax, delta |
|---|
| 37 | integer j |
|---|
| 38 | real(8) f, g, gprime |
|---|
| 39 | |
|---|
| 40 | if (nz_max<nz_delta .or. nz_delta<=1) then |
|---|
| 41 | stop 'inappropriate input to smartzfac' |
|---|
| 42 | endif |
|---|
| 43 | f = 1.05 |
|---|
| 44 | do j=1,7 ! Newton iteration |
|---|
| 45 | !print *,j,f |
|---|
| 46 | g = (f-3+2*f**(nz_max-1))/(f-3+2*f**(nz_delta-1))-zmax/delta |
|---|
| 47 | gprime=(1+2*(nz_max-1)*f**(nz_max-2))/(f-3+2*f**(nz_delta-1)) - & |
|---|
| 48 | & (f-3+2*f**(nz_max-1))*(1+2*(nz_delta-1)*f**(nz_delta-2))/ & |
|---|
| 49 | & (f-3+2*f**(nz_delta-1))**2 |
|---|
| 50 | f = f-g/gprime |
|---|
| 51 | enddo |
|---|
| 52 | smartzfac = f |
|---|
| 53 | if (smartzfac<1. .or. smartzfac>2.) then |
|---|
| 54 | print *,'zfac=',smartzfac |
|---|
| 55 | stop 'unwanted result in smartzfac' |
|---|
| 56 | endif |
|---|
| 57 | end function smartzfac |
|---|
| 58 | |
|---|
| 59 | |
|---|
| 60 | |
|---|
| 61 | !-grid-dependent utility functions |
|---|
| 62 | |
|---|
| 63 | function colint(y,z,nz,i1,i2) |
|---|
| 64 | ! column integrates y |
|---|
| 65 | implicit none |
|---|
| 66 | integer, intent(IN) :: nz, i1, i2 |
|---|
| 67 | real(8), intent(IN) :: y(nz), z(nz) |
|---|
| 68 | real(8) colint |
|---|
| 69 | integer i |
|---|
| 70 | real(8) dz(nz) |
|---|
| 71 | |
|---|
| 72 | dz(1) = (z(2)-0.)/2 |
|---|
| 73 | do i=2,nz-1 |
|---|
| 74 | dz(i) = (z(i+1)-z(i-1))/2. |
|---|
| 75 | enddo |
|---|
| 76 | dz(nz) = z(nz)-z(nz-1) |
|---|
| 77 | colint = sum(y(i1:i2)*dz(i1:i2)) |
|---|
| 78 | end function colint |
|---|
| 79 | |
|---|
| 80 | |
|---|
| 81 | |
|---|
| 82 | subroutine heatflux_from_temperature(nz,z,T,k,H) |
|---|
| 83 | ! calculates heat flux from temperature profile |
|---|
| 84 | ! like k, the heat flux H is defined mid-point |
|---|
| 85 | implicit none |
|---|
| 86 | integer, intent(IN) :: nz |
|---|
| 87 | real(8), intent(IN) :: z(nz), T(nz), k(nz) |
|---|
| 88 | real(8), intent(OUT) :: H(nz) |
|---|
| 89 | integer j |
|---|
| 90 | |
|---|
| 91 | ! H(1) = -k(1)*(T(1)-Tsurf)/z(1) |
|---|
| 92 | H(1) = 0. ! to avoid ill-defined value |
|---|
| 93 | do j=2,nz |
|---|
| 94 | H(j) = -k(j)*(T(j)-T(j-1))/(z(j)-z(j-1)) |
|---|
| 95 | enddo |
|---|
| 96 | end subroutine heatflux_from_temperature |
|---|
