Listing 3: FORTRAN Program fdtrns.for, Finite Difference Method to Solve [24] and [32] Note: Subroutines tridia and wprop omitted c fdtrns implicit double precision (a-h, k, o-z) parameter (iz=500) dimension aa(0:iz), bb(0:iz), cc(0:iz), dd(0:iz) dimension th(0:iz), h(0:iz), ht(0:iz), k(0:iz), z(0:iz) dimension cl(0:iz), u(0:iz), v(0:iz) open (1, file='outtrnsd.dat', status='unknown') zl = .1d0 tl = 36000.d0 dt = 10.d0 np = 100 ae = 0.5d0 de = 2.d-10 co = 0.1d0 dtp = 3600.d0 tp = dtp dz = zl/dble(np) rz = dt/dz rzz = dt/(dz*dz) write (1, *) np write (1, '(f10.4)') zl, ae, dz, tl, dt, dtp, de, co c initialize to hydrostatic conditions do i = 0, np z(i) = -zl + dble(i)*dz ht(i) = -zl h(i) = -dble(i)*dz p = -h(i) call wprop (p, th(i), c, k(i)) end do knp1 = k(np) write (*, *) c do Picard iterations until solution or max it do it = 1, 100 do i = 1, np-1 km = 0.5d0*(k(i-1)+k(i)) kp = 0.5d0*(k(i)+k(i+1)) aa(i) = -km bb(i) = km+kp cc(i) = -kp dd(i) = 0.d0 end do dd(1) = -0.5d0*(k(0)+k(1))*zl km = 0.5d0*(k(np-1)+k(np)) hnp1 = h(np-1) - 2.d0*dz*(1.d0+ae) p = -hnp1 call wprop (p, thnp1, cnp1, knp1) kp = 0.5d0*(k(np)+knp1) aa(np) = -(km+kp) bb(np) = km + kp cc(np) = 0.d0 dd(np) = -2.d0*dz*ae*kp call tridia (np, aa, bb, cc, dd) do i = 1, np ht(i) = dd(i) h(i) = ht(i) - z(i) p = -h(i) call wprop (p, th(i), c, k(i)) end do c quit if maximum relative difference in cell flows is < 1.e-6 qmin = -0.5d0*(k(1)+k(0))*(ht(1)-ht(0))/dz qmax = qmin do i = 2, np q = -0.5d0*(k(i-1)+k(i))*(ht(i)-ht(i-1))/dz if (q .lt. qmin) qmin = q if (q .gt. qmax) qmax = q end do rf = (qmax-qmin)/(ae*k(np)) write (*, '(1h, i5, f13.6)') it, rf if (rf .lt. 1.d-6) go to 100 end do 100 continue write (*, *) qmin, k(np), qmin/k(np) write (*, *) c write out profile data write (1, *) write (1, '(g11.4)') qmin, qmax, rf write (1, *) do i = 0, np write (1, '(f7.4,f9.3,f7.4,g11.4)') z(i), h(i), th(i), k(i) end do write (*, *) qe = 0.5d0*(qmin+qmax) do i = 0, np-1 c v(i) = qe/(0.5d0*(th(i+1)+th(i))) v(i) = qe/th(i) end do do i = 0, np cl(i) = co u(i) = th(i)*co end do t = 0.d0 do while ((t+dt) .le. tl) t = t + dt write (*, '(1h+, f10.2)') t do i = 1, np-1 tdm = 0.5d0*(th(i-1)+th(i))*de tdp = 0.5d0*(th(i)+th(i+1))*de aa(i) = -(rz*v(i-1) + rzz*tdm/th(i-1)) bb(i) = 1.d0 + rz*v(i) + rzz*(tdm+tdp)/th(i) cc(i) = -rzz*tdp/th(i+1) dd(i) = u(i) end do dd(1) = dd(1) - aa(1)*u(0) tdm = tdp aa(np) = -(rz*v(np-1) + rzz*tdm/th(np-1)) bb(np) = 1.d0 + rzz*tdm/th(np) cc(np) = 0.d0 dd(np) = u(np) call tridia (np, aa, bb, cc, dd) do i = 1, np u(i) = dd(i) cl(i) = u(i)/th(i) end do if (t .ge. tp) then write (1, *) do i = 0, np write (1, '(f7.4,2g11.4)') z(i), u(i), cl(i) end do tp = tp + dtp end if end do close (1) end