Module ETdata
   implicit none
   save
   character(len=80) :: fn
   integer :: year, month, nd, ios, j, im, nm, nn(12)
   real :: phid, h, p, pa, pps, azd, td, t, s, v, vd, delta, &
theta, tc, c0, c1, c2, w, phi, rh
real :: alpha(2) = (/17.27, 21.88/), beta(2) = (/237.3, 265.5/),&
pi = 3.14159254, gamma(2), ftz(2), fz(2), az, azz, cz, &
			cosz, comega, eta, ge, z, tmp, a0, g0, b, rt, rtc, &
			zeta, ft, lambda, et, etp, tau, rho, delp, delt, vp, &
tp, es, d, salt, omega, taua, g, a, b0, rtp, etw, b1, &
b2, ltheat, xm(4,12), gdew, rw, ep, ew
end Module ETdata

program ET
!  Calculates Morton's ET
!  05 April 2012
   use ETdata
   implicit none
   integer :: flag
   open (10, file = 'ETpaper.par', status = 'old')
   read (10, 100) fn; 100 format(a)
	open (11, file = fn, status = 'old')										
! 	Input data file
	read (10, 100) fn
	open (21, file = fn)													
! 	Output file
	read (11, 100) fn
	write (21, 100) fn
   read (10, *) flag
   if (flag == 1) then
!		b0 = 1.; b1 = 14.; b2 = 1.20; fz(1) = 28.	 	! Morton's values
		b0 = 1.; b1 = 13.4; b2 = 1.13; fz(1) = 29.2		! QJ's values
		es = 5.22e-8
		write (21, '("fz, b1, and b2", 3f8.2)'), fz(1), b1, b2
		read (11, *) phid, h, pa, nm
		write (21, 200) phid, h, pa, nm
		200 format('Latitude =', f7.2, '  Station height =', f8.2, &
' m  Mean annual precipitation =', f8.1, ' mm'/  &
'#  of months =', i5 )
		write (21, 201)
		201 format(' Year Month  Temp     Rel Hum  Sunshine'  &
'Net Rad     PET        WET        AET     ft'/ &
			'             (oC)       (%)    (Ratio)   (mm)' &
'       (mm)       (mm)       (mm)')
	else
		b0 = 1.12
		b1 = 13.
		b2 = 1.12
		es = 5.5e-8
		fz(1) = 25.
		if (flag == 2) then
			read (11, *) phid, h, pa, nm
			write (21, 200) phid, h, pa, nm
		else if (flag == 3 ) then
			read (11, *) phid, h, pa, nm, d, salt
			write (21, 202) phid, h, pa, d, salt, nm
			202 format('Latitude =', f7.2, '  Station height =', &
 f8.2, ' m', '  Mean annual precipitation =', &
				 f8.1, ' mm'/ ' Lake depth =', f7.2, ' m',    &
 '  Salt =', f8.1, ' ppm', ' #  of months =', i5 )
		else if (flag == 4 ) then
			read (11, *) phid, h, pa, nm, salt
			write (21, 203) phid, h, pa, salt, nm
			203 format('Latitude =', f7.2, '  Station height =', &
f8.2, ' m', '  Mean annual precipitation =', &
				f8.1, ' mm'/ '  Salt =', f8.1, ' ppm',       &
' #  of months =', i5 )
		end if
	end if
	read (11, *)														
! 	Skip a line
	phi = phid*pi/180.
!  	Compute the ratio of atmospheric pressure at the station to that
!  at sea level
	pps = ((288. - 0.0065*h)/288)**5.256							  !C-1
	p = pps*1013.
!  	Estimate the zenith value of the dry season snow free clear sky
!  albedo
azd = 0.26 - 0.00012*pa*(1 + abs(phid/42) +       &
 (phid/42)**2)*pps**0.5									  !C-2
	if (azd < 0.11) azd = 0.11
	if (azd > 0.17) azd = 0.17
	gamma(1) = 0.66*pps
	gamma(2) = gamma(1)/1.15
	fz(2) = fz(1)*1.15
	ftz(1) = fz(1)*sqrt(1./pps)
	ftz(2) = ftz(1)*1.15
	nn = 0 ; xm =0.
	if (flag == 1) then
		call crae(1)
	else if (flag == 2) then
		write (21, 204)
		204 format(' Year Month  Temp     Rel Hum  Sunshine  ',  &
 'Net Rad     PET        AET'/ &
			 '             (oC)       (%)    (Ratio)   (mm)',  &
  '       (mm)       (mm)')
		call crwe(2)
	else if (flag == 3) then
		write (21, 204)
		call crle1(3)
	else if (flag == 4) then
		write (21, 205)
		205 format(' Year Month  Temp     Rel Hum  Sunshine',  &
 '    Depth   Net Rad     PET        AET'/ &
					 '             (oC)       (%)    (Ratio)', &
 '      (m)     (mm)       (mm)       (mm)')
		call crle2(4)
	end if
end program ET

subroutine crae(flag)
!  	Calculates the actual evapotranspiration from landscape
   	use ETdata
   	implicit none
	integer :: flag
	do im = 1, nm
		read (11, *, iostat = ios) year, month, nd, t, rh, s
		if (ios < 0) exit
		if (s < 0. .OR. rh < 0. .OR. t < -999.) cycle
		call sub1(flag)
!		Estimate the net radiation for soil-plant surfaces at air
!		temperature (rt), the stability factor (sf), the vapour
!		pressure coefficient (ft) and the heat transfer 
!		coefficient (lambda)
		rt = (1. - a)*g - b											!C-27
		rtc = rt														!C-28
		if (rtc < 0.) rtc = 0.										!C-28a
		zeta = 1./(0.28*(1. + vd/v) + delta*rtc*pps**0.5/ &
(gamma(j)*b0*fz(j)*(v - vd)))						!C-29
		if (zeta < 1.) zeta = 1.										!C-29a
		ft = fz(j)/zeta/pps**0.5										!C-30
		lambda = gamma(j) + 4.*es*(t + 273)**3/ft					!C-31
!		Estimate the final values from the following quickly
!		converging iterative solution of the vapour transfer
!		and energy balance equations
		tp = t															
		vp = v
		delp = delta
		do
			delt = (rt/ft + vd - vp +lambda*(t - tp))/  &
(delp + lambda)								!C-32
			tp = tp + delt											!C-33
			vp = 6.11*exp(alpha(j)*tp/(tp + beta(j)))				!C-34
			delp = alpha(j)*beta(j)*vp/(tp + beta(j))**2			!C-35
			if (abs(delt) < 0.01) exit
		end do
		
!		Estimate the potential evapotranspiration (etp), the net
!		radiation for soil-plant surfaces at equilibrium temperature
!		(rtp) and the wet-environment areal evapotranspiration (etw)
		etp = rt - lambda*ft*(tp - t)								!C-36
		rtp = etp + gamma(j)*ft*(tp - t)							!C-37
		etw = b1 + b2*rtp/(1. + gamma(j)/delp)					!C-38
		if (etw < etp/2.) etw = etp/2.								!C-38a
		if (etw > etp) etw = etp										!C-38a
!		Estimate the areal evapotranspiration (et) from the
!		complementary relationship
		et = 2.*etw - etp												!C-39
		ltheat = 28.5
		if (t < 0.) ltheat = 28.5*1.15
		et = nd*et/ltheat
		rt = nd*rt/ltheat
		etp = nd*etp/ltheat
		etw = nd*etw/ltheat
		write (21, 201) year, month, t, rh, s, rt, etp, etw, et, ft
		201 format(i5, i3, 8f10.2)
		xm(1,month) = xm(1,month) + rt
		xm(2,month) = xm(2,month) + etp
		xm(3,month) = xm(3,month) + etw
		xm(4,month) = xm(4,month) + et
		nn(month) = nn(month) + 1
	end do
	write (21, 202)
	202 format(/10x,'Mean Values'/'Month   Net Rad      PET'  & 
'       WET       AET  # of months')
	do month = 1, 12
		xm(:,month) = xm(:,month)/nn(month)
		write (21, 203) month, (xm(j,month), j = 1, 4), nn(month)
		203 format(i5, 4f10.2, i5)
	end do
	write (21, 204) (sum(xm(j,:)), j = 1, 4)
	204 format('Annual', f9.1, 3f10.2)
end subroutine crae

subroutine sub1(flag)
	use ETdata
	implicit none
	integer :: flag
	if (t < 0.) then
		j = 2
	else
		j = 1
	end if
!	Compute the saturation vapour pressure at td (vd), at t (v) and
!	the slope of the saturation vapour pressure at t (delta)
	v = 6.11*exp(alpha(j)*t/(t + beta(j)))						 !C-4
!	Calculate dew point temperature from relative humidity 
!	(Lawrence 2005)
	gdew = 17.625*t/(243.04 + t) + log(rh/100.)
	td = 243.04*gdew/(17.625 - gdew)
	vd = 6.11*exp(17.27*td/(td + 237.3))
	delta = alpha(j)*beta(j)*v/(t + beta(j))**2					 !C-5
!	Compute various angles and functions leading up to an estimate 
!	of the extra-atmospheric global radiation (ge)
	theta = 23.2*sin((29.5*month - 94.)*pi/180.)*pi/180.		 !C-6
	z = phi - theta
	cosz = cos(z)													  	 !C-7
	if (cosz < 0.001) cosz = 0.001									 !C-7a
	if (cosz > 1.) cosz = 1.
	z = acos(cosz)
	comega = 1. - cosz/cos(phi)/cos(theta)	 					 !C-8
	if (comega < -1.) comega = -1.									 !C-8a
	omega = acos(comega)
	cz = cosz + (sin(omega)/omega - 1.)*cos(phi)*cos(theta)		 !C-9
	eta = 1. + sin((29.5*month - 106)*pi/180.)/60.				 !C-10
	ge = 1354*omega*cz/pi/eta/eta									!C-11
!	Estimate the zenith value of snow-free clear-sky albedo (azz),
!	zenith value of clear-sky albedo (az) and the clear sky 
!	albedo (a0)
	if (flag > 1) then
		azz = 0.05
	else
		azz = azd														!C-12
		if (azz > 0.5*(0.91 - vd/v)) azz = 0.5*(0.91 - vd/v)	!C-12A
	end if					
	c0 = v - vd														!C-13
	if (c0 < 0.) c0 = 0.												!C-13a
	if (c0 > 1.) c0 = 1.												!C-13a
	az = azz + (1. - c0*c0)*(0.34 - azz)							!C-14
	a0 = az*(exp(1.08) - (2.16*cosz/pi +   &
sin(z))*exp(0.012*z*180./pi))/1.473/(1. - sin(z))			!C-15
!	Estimate precipitable water (w) and a turbidity coefficient (tc)
	w = vd/(0.49 + t/129.)											!C-16
	c1 = 21 - t														!C-17
	if (c1 < 0.) c1 = 0.												!C-17a
	if (c1 > 5.) c1 = 5.												!C-17a
	tc  = (0.5 + 2.5*cz*cz)*exp(c1*(pps - 1.))					!C-18
!	Compute the transmittancy of clear skies to direct beam 
!	solar radiation (tau)
	tau = exp(-0.089*(pps/cz)**0.75 - 0.083*(tc/cz)**0.90 - &
0.029*(w/cz)**0.60)	     								!C-19
!	Estimate the part of tau that is the result of absorption (taua)
	taua = exp(-0.0415*(tc/cz)**0.90 - (0.0029**0.5)*  &
(w/cz)**0.3)											!C-20
	tmp = exp(-0.0415*(tc/cz)**0.90 - 0.029*(w/cz)**0.6)
	if (taua < tmp) taua = tmp										!C-20a
!	Compute the clear-sky global radiation (g0) and the incident
!	global radiation (g)
	g0 = ge*tau*(1. + (1. - tau/taua)*(1. + a0*tau))	 			!C-21
	s = s*pi/24./omega			
	if (s > 1.) s = 0.999
	g = s*g0 + (0.08 + 0.30*s)*(1. - s)*ge						!C-22
!	Estimate the average albedo (a)
	a = a0*(s + (1. - s)*(1. - z*180/pi/330))						!C-23
!	Estimate the proportional increase in atmospheric radiation
!	due to clouds (rho)
	c2 = 10.*(vd/v - s - 0.42)										!C-24
	if (c2 < 0.) c2 = 0.												!C-24a
	if (c2 > 1.) c2 = 1.												!C-24a
	rho = 0.18*((1. - c2)*(1. - s)**2 + c2*(1. - s)**0.5)/pps	!C-25
!	Calculate the net long-wave radiation loss for soil-plant
!	surface at air temperature (b)
	b = es*(t + 273)**4*(1. - (0.71 + 0.007*vd*pps)*(1. + rho)) !C-26
	tmp = 0.05*es*(t + 273)**4
	if ( b < tmp) b = tmp											!C-26a
end subroutine sub1

subroutine crwe(flag)
!  	Calculates Morton's ET - Shallow lake evaporation
!	05 April 2012
	use ETdata
	implicit none
	integer :: flag
	do im = 1, nm
		read (11, *, iostat = ios) year, month, nd, t, rh, s
		if (ios < 0) exit
		if (s < 0. .OR. rh < 0. .OR. t < -999.) cycle
		call sub1(flag)
!		Estimate the net radiation for water surface at air
!		temperature (rw),
!		the stability factor (sf), the vapour pressure 
!		coefficient (ft) and the heat transfer coefficient (lambda)
		rw = (1. - a)*g - b
		rtc = rw
		if (rtc < 0.) rtc = 0.
		zeta = 1./(0.28*(1. + vd/v) + delta*rtc*pps**0.5/  &
(gamma(j)*b0*fz(j)*(v - vd)))
		if (zeta < 1.) zeta = 1.
		ft = fz(j)/zeta/pps**0.5
		lambda = gamma(j) + 4.*es*(t + 273)**3/ft
!		Estimate the final values from the following quickly 
!		converging iterative solution of the vapour transfer
!		and energy balance equations
		tp = t
		vp = v
		delp = delta
		do
			delt = (rw/ft + vd - vp +lambda*(t - tp))/(delp + lambda)
			tp = tp + delt
			vp = 6.11*exp(alpha(j)*tp/(tp + beta(j)))
			delp = alpha(j)*beta(j)*vp/(tp + beta(j))**2
			if (abs(delt) < 0.01) exit
		end do
		
!		Estimate the potential evaporation (ep), the net radiation
!		for soil-plant surfaces at equilibrium temperature (rtp)
!		and the shallow lake evaporation (ew)
		ep = rw - lambda*ft*(tp - t)
		rtp = ep + gamma(j)*ft*(tp - t)
		ew = b1 + b2*rtp/(1. + gamma(j)/delp)
		if (ew < ep/2.) ew = ep/2.
		if (ew > ep) ew = ep
		ltheat = 28.5
		if (t < 0.) ltheat = 28.5*1.15
		ep = nd*ep/ltheat
		rw = nd*rw/ltheat
		ew = nd*ew/ltheat
		write (21, 201) year, month, t, rh, s, rw, ep, ew
		201 format(i5, i3, 6f10.2)
		xm(1,month) = xm(1,month) + rw
		xm(2,month) = xm(2,month) + ep
		xm(3,month) = xm(3,month) + ew
		nn(month) = nn(month) + 1
	end do
	write (21, *) ' Mean values'
	write (21, *) 'Month  Net Rad     PET        AET'
	do month = 1, 12
		xm(:,month) = xm(:,month)/nn(month)
		write (21, 203) month, (xm(j,month), j = 1, 3)
		203 format(i5, 4f10.2)
	end do
	write (21, 204) (sum(xm(j,:)), j = 1, 3)
	204 format('Annual', f9.1, 3f10.2)
end subroutine crwe

subroutine crle1(flag)
!  	Calculates Morton's ET - Deep lake evaporation
!	05 April 2012
	use ETdata
	implicit none
	integer, parameter :: nmx = 1212
	integer :: ic, iy(nmx), mon(nmx), tnd(nmx), flag
	real ::  tgw(nmx), tv(nmx), tvd(nmx), ts(nmx), gl(nmx), &
tt(nmx), trh(nmx), tdelta(nmx)
	do im = 1, nm
		read (11, *, iostat = ios) year, month, nd, t, rh, s
		if (ios < 0) exit
		if (s < 0. .OR. rh < 0. .OR. t < -999.) cycle
		call sub1(flag)
		tgw(im+12) = (1. - a)*g - b
		tt(im) = t
		trh(im) = rh
		ts(im) = s
		tv(im) = v
		tvd(im) = vd
		tdelta(im) = delta
		iy(im) = year
		tnd(im) = nd
		mon(im) = month
	end do
	call deepLake(nmx, nm, tgw, gl, d, salt)
	do im = 1, nm
		t = tt(im)
		rh = trh(im)
		s = ts(im)
		v = tv(im)
		vd = tvd(im)
		delta = tdelta(im)
! 	   	Estimate the net radiation for water surface at air
!		temperature (rw),
!		the stability factor (sf), the vapour pressure coefficient
!		(ft) and the heat transfer coefficient (lambda)
		rw = gl(im)
		rtc = rw
		if (rtc < 0.) rtc = 0.
		zeta = 1./(0.28*(1. + vd/v) + delta*rtc*pps**0.5/   &
(gamma(j)*b0*fz(j)*(v - vd)))						!C-29
		if (zeta < 1.) zeta = 1.										!C-29a
		ft = fz(j)/zeta/pps**0.5										!C-30
		lambda = gamma(j) + 4.*es*(t + 273)**3/ft					!C-31
!	   	Estimate the final values from the following quickly
!		converging iterative solution of the vapour transfer
!		and energy balance equations
		tp = t
		vp = v
		delp = delta
		ic = 0
		do
			delt = (rw/ft + vd - vp +lambda*(t - tp))/   & 
(delp + lambda) 															!C-32
			tp = tp + delt													!C-33
			vp = 6.11*exp(alpha(j)*tp/(tp + beta(j)))				!C-34
			delp = alpha(j)*beta(j)*vp/(tp + beta(j))**2			!C-35
			if (abs(delt) < 0.01) exit												
			ic = ic + 1
			if (ic > 100) then
				print *, 'Did not converge for month', im
				exit
			end if
		end do
!
!	   	Estimate the potential evaporation (ep), the net radiation 
!		for soil-plant surfaces at equilibrium temperature (rtp) and
!		the deep lake evaporation (ew)
		ep = rw - lambda*ft*(tp - t)
		rtp = ep + gamma(j)*ft*(tp - t)
		ew = b1 + b2*rtp/(1. + gamma(j)/delp)
		if (ew < ep/2.) ew = ep/2.
		if (ew > ep) ew = ep
		ltheat = 28.5
		if (t < 0.) ltheat = 28.5*1.15
		nd = tnd(im)
		ep = nd*ep/ltheat
		rw = nd*rw/ltheat
		ew = nd*ew/ltheat
		write (21, 201) iy(im), mon(im), t, rh, s, rw, ep, ew
		201 format(i5, i3, 6f10.2)
		xm(1,mon(im)) = xm(1,mon(im)) + rw
		xm(2,mon(im)) = xm(2,mon(im)) + ep
		xm(3,mon(im)) = xm(3,mon(im)) + ew
		nn(mon(im)) = nn(mon(im)) + 1
	end do
	write (21, *) ' Mean values'
	write (21, *) 'Month  Net Rad     PET        AET'
	do month = 1, 12
		xm(:,month) = xm(:,month)/nn(month)
		write (21, 203) month, (xm(j,month), j = 1, 3)
		203 format(i5, 4f10.2)
	end do
	write (21, 204) (sum(xm(j,:)), j = 1, 3)
	204 format('Annual', f9.1, 3f10.2)
end subroutine crle1

subroutine deepLake(nmx, nm, tgw, gl, d, salt)
!  	Calculates the available solar and waterborne heat(gl) using
!	Morton 1986
	implicit none
	integer :: im, nmx, nm, ti, t1, i, i1, ii
	real :: gw(nmx), tgw(nmx), gl(nmx), t0, t, glb, gle, &
salt, d, f, k
	t0 = 0.96 + 0.013*d												!(6)
	if (t0 < 0.039*d) t0 = 0.039*d
	if (t0 > 0.13*d) t0 = 0.13*d
	t = t0/(1. + salt/27000)**2											!(7)
	if (t > 6.) t = 6.
	k = t0/(1. + (d/93)**7)											!(8)
!  	Calculates the delayed solar and waterborne energy (gw)
	ti = int(t)
	t1 = ti + 1
	f = t - ti
	i = 12
	tgw(1:12) = tgw(13:24)
	do im = 1, nm
		i = i + 1
		i1 = i - t1
		ii = i - ti
		gw(im) = tgw(ii) + f*(tgw(i1) - tgw(ii))							!(3)
	end do
! 	Calculates the available solar and waterborne heat(gl)
	glb = 50.
	do i = 1, 2	
		do im = 1, 12
			gle = glb + (gw(im) - glb)/(k + 0.5)					!(4)
			gl(im) = (glb + gle)/2.									!(5)
			glb = gle
		end do
	end do
	do im = 1, nm
		gle = glb + (gw(im) - glb)/(k + 0.5)								!(4)
		gl(im) = (glb + gle)/2.										!(5)
		glb = gle
	end do
end subroutine deepLake

subroutine crle2(flag)
!  	Calculates Morton's ET - Deep lake evaporation
!	05 April 2012
	use ETdata
	implicit none
	integer, parameter :: nmx = 1212
	integer :: ic, iy(nmx), mon(nmx), tnd(nmx), flag
	real ::  tgw(nmx), tv(nmx), tvd(nmx), ts(nmx), gl(nmx), &
tt(nmx), trh(nmx), tdelta(nmx), da(nmx)
	do im = 1, nm
		read (11, *, iostat = ios) year, month, nd, t, rh, s, da(im)
		if (ios < 0) exit
		if (s < 0. .OR. rh < 0. .OR. t < -999.) cycle
		call sub1(flag)
		tgw(im+12) = (1. - a)*g - b
		tt(im) = t
		trh(im) = rh
		ts(im) = s
		tv(im) = v
		tvd(im) = vd
		tdelta(im) = delta
		iy(im) = year
		mon(im) = month
		tnd(im) = nd
	end do
	call deepLakeVD(nmx, nm, tgw, gl, da, salt)
	do im = 1, nm
		t = tt(im)
		rh = trh(im)
		s = ts(im)
		v = tv(im)
		vd = tvd(im)
		delta = tdelta(im)
! 	  	Estimate the net radiation for water surface at air
!		temperature (rw),
!		the stability factor (sf), the vapour pressure coefficient 
!		(ft) and the heat transfer coefficient (lambda)
		rw = gl(im)
		rtc = rw
		if (rtc < 0.) rtc = 0.
		zeta = 1./(0.28*(1. + vd/v) + delta*rtc*pps**0.5/  &
 (gamma(j)*b0*fz(j)*(v - vd)))						!C-29
		if (zeta < 1.) zeta = 1.										!C-29a
		ft = fz(j)/zeta/pps**0.5										!C-30
		lambda = gamma(j) + 4.*es*(t + 273)**3/ft					!C-31
!	   	Estimate the final values from the following quickly
! 		converging iterative solution of the vapour transfer and 
!		energy balance equations
		tp = t
		vp = v
		delp = delta
		ic = 0
		do
			delt = (rw/ft + vd - vp +lambda*(t - tp))/(delp + lambda)			!C-32
			tp = tp + delt											!C-33
			vp = 6.11*exp(alpha(j)*tp/(tp + beta(j)))				!C-34
			delp = alpha(j)*beta(j)*vp/(tp + beta(j))**2			!C-35
			if (abs(delt) < 0.001) exit											! ****************
			ic = ic + 1
			if (ic > 100) then
				print *, 'Did not converge for month', im
				exit
			end if
		end do
	
!	   	Estimate the potential evaporation (ep), the net radiation
!		for soil-plant
!		surfaces at equilibrium temperature (rtp) and the deep lake
!		evaporation (ew)
!
		ep = rw - lambda*ft*(tp - t)
		rtp = ep + gamma(j)*ft*(tp - t)
		ew = b1 + b2*rtp/(1. + gamma(j)/delp)
		if (ew < ep/2.) ew = ep/2.
		if (ew > ep) ew = ep
		ltheat = 28.5
		if (t < 0.) ltheat = 28.5*1.15
		nd = tnd(im)
		ep = nd*ep/ltheat
		rw = nd*rw/ltheat
		ew = nd*ew/ltheat
		write (21, 201) iy(im), mon(im), t, rh, s, da(im), rw, ep, ew
		201 format(i5, i3, 7f10.2)
		xm(1,mon(im)) = xm(1,mon(im)) + rw
		xm(2,mon(im)) = xm(2,mon(im)) + ep
		xm(3,mon(im)) = xm(3,mon(im)) + ew
		nn(mon(im)) = nn(mon(im)) + 1
	end do
	write (21, *) ' Mean values'
	write (21, *) 'Month  Net Rad     PET        AET'
	do month = 1, 12
		xm(:,month) = xm(:,month)/nn(month)
		write (21, 203) month, (xm(j,month), j = 1, 3)
		203 format(i5, 4f10.2)
	end do
	write (21, 204) (sum(xm(j,:)), j = 1, 3)
	204 format('Annual', f9.1, 3f10.2)
end subroutine crle2
		
subroutine deepLakeVD(nmx, nm, tgw, gl, da, salt)
!  	Calculates the available solar and waterborne heat(gl) using
!	Morton 1986
	implicit none
	integer :: im, nmx, nm, ti, t1, i, i1, ii
	real :: gw(nmx), tgw(nmx), gl(nmx), t0, t, glb, gle, salt, &
d, f, da(nmx), ka(nmx)
	tgw(1:12) = tgw(13:24)
	i = 12
	do im = 1, nm
		d = da(im)
		t0 = 0.96 + 0.013*d											!(6)
		if (t0 < 0.039*d) t0 = 0.039*d
		if (t0 > 0.13*d) t0 = 0.13*d
		t = t0/(1. + salt/27000)**2									!(7)
		if (t > 6.) t = 6.
		ka(im) = t0/(1. + (d/93)**7)										!(8)
!  	Calculates the delayed solar and waterborne energy (gw)
		ti = int(t)
		t1 = ti + 1
		f = t - ti
		i = i + 1
		i1 = i - t1
		ii = i - ti
		gw(im) = tgw(ii) + f*(tgw(i1) - tgw(ii))					!(3)
	end do
!  	Calculates the available solar and waterborne heat(gl)
	glb = 50.
	do i = 1, 2	
		do im = 1, 12
			gle = glb + (gw(im) - glb)/(ka(im) + 0.5)					!(4)
			gl(im) = (glb + gle)/2.									!(5)
			glb = gle
		end do
	end do
	do im = 1, nm
		gle = glb + (gw(im) - glb)/(ka(im) + 0.5)						!(4)
		gl(im) = (glb + gle)/2.										!(5)
		glb = gle
	end do
end subroutine deepLakeVD


