subroutine moblam8
!=========================
!....version 10-feb-99, Peter Schluessel
!...surface renewal model including solar and longwave radiative
!...heat fluxes as well as turbulent fluxes of sensible and latent heat
!...and impact of rain on the cool skin and air-sea gas exchange
!...algorithm descriptions in:
!...Soloviev, A.V., P. Schluessel (1994): JPO, 24, 1339-1346 + 1965(errata).
!...Soloviev, A.V., P. Schluessel (1996): BLM, 77, 45-68.
!...Schluessel, P., A.V. Soloviev, W.J. Emery (1997): BLM, 82, 437-472.
!...Craeye, C., P. Schluessel (1998): BLM, 89,349-355.
!
!...If comments and program code disagree probably both are wrong.
!
!...Alexander Soloviev, PPSIO, Matlab version
!...Peter Schluessel, Met. Inst., University Hamburg        18-jan-94
!...last change of program code before public release	     8-nov-94
!...inclusion of surface and volume heat flux due to rain   20-feb-96
!...inclusion of fresh-water skin                           20-feb-96
!...inclusion of friction velocity related to kinetic energy
!...of rain                                                 29-feb-96
!...inclusion of toga-coare bulk flux algorithm             29-feb-96
!...horizontal surface stress produced by rainfall           7-mar-96
!...adapted to fortran90        		                    20-mar-97
!...inclusion of renewal bulk flux algorithm		        11-jun-97
!...minor changes                                           26-jun-97
!...major changes in txrain and rfric                       26-mar-98
!...inclusion of bubble-mediated gas-transfer velocity      19-oct-98
!...modification of renewal model for the air               20-oct-98
!...minor modifications					    10-feb-99 
!
!...externals: skin_ln, t_star, rain_ln, frain_ln, txrain, erf, gamma_p,
!              hmg, gseries, gamma_ln, gamma_inc, gamma_q, psyceq, densit,
!              spec_h, viscos, cdepth, rfric, gammain, ren_bulk
!              humidity, rstress, rain_z0, c_ustar, bubble_k
!
!...compile with double precision floating variables		!!!!!!!!!!!!
!
!-----------------------------------------------------------------------
!
!...input: qlong, longwave net flux (W/m^2)
!          qsol, shortwave net flux (W/m^2)
!          ws, wind speed (m/s)
!          tsurf, surface temperature (K)
!          sal, salinity in parts per thousand
!          p, surface air pressure (hPa)
!          t, surface air temperature (K)
!          q, specific humidity (kg/kg)
!          rr, surface rainrate (mm/h)
!	   zu, anemometer height (m)
!	   zt, thermometer height (m)
!	   zq, hygrometer height (m)
!          skin_flag, true if tsurf is skin temperature, false if it is bulk temperature
!...output: tmax, renewal time in water (s)
!           ustar_a, friction velocity in air (m/s)
!           ustar_w, friction velocity due to wind stress (m/s)
!           ustar_h, friction velocity by horizontal rain sress (m/s)
!           ustar_r, friction velocity due to rain-induced mixing (m/s)
!           ustar, total friction velocity in water(m/s)
!           tr, rain temperature (K)
!           tc, skin cooling/warming by surface energy fluxes (K)
!           ts, skin warming by solar radiation (K)
!           trs, skin cooling by rain cumulating at surface (K)
!           trv, skin cooling by volume flux due to rain (K)
!           srs, fresh water skin by rain cumulating at surface (parts per thousand)
!           srv, fresh water skin by volume flux (parts per thousand)
!			gc, gas-transfer velocity for CO2 (m/s)
!			bc, bubble-mediated gas-transfer velocity for CO2 (m/s) 
!
 	save
 	real, parameter :: t0=273.15, rgas=287.05
	real ws,qq,p
	real tau,ustar,qstar,t_star
	real rr,zl,zo
	real zot,zoq
	logical skin_flag
	external rstress,densit,t_star,frain_ln,rain_ln,skin_ln
	common/skin_in/ws,t,q,tsurf,qs,qsol,qlong,rain,sal,zu,zt,zq,p,skin_flag
    common/skin_out/tmax,ustar_w,ustar_h,ustar_r,ustar,tr,tc,ts,trs,trv,srs,srv,gc,bc
	common/bulk_out/ustar_a,tstar,qstar,zo,zot,zoq,zl, tau, hf, ef, tren
	data zu/10./,zt/10./,zq/10./

!...add cheat value to rainrate to avoid numerical instabilities at rain=0.
 	rain = rain + 1.d-6
	rr = rain

! call to bulk-algorithm 
	call ren_bulk
 
! calculate dewpoint temperature
 	ta=t
	ttd=dewpt(q,p)
	td=min(ttd-t0,t)

!...water-vapour pressure in hpa
 	ee=goff_gratch(td+t0)

!...specific humidity in g/kg
 	qq=622.d0*ee/(p-0.378d0*ee)
!
!...heat flux in upper ocean
	q0=qlong+ef+hf
	ql=ef

!...air density
 	roha=p*100.d0/((ta+t0)*rgas)

!...ocean water density
 	rohw=densit(tsurf,sal)

!...wind-induced friction velocity in upper ocean
 	ustar_w=ustar_a*sqrt(roha/rohw)

!...surface stress produced by rainfall
 	tau_h=rstress(ws,rr)

!...friction velocity related to it
 	ustar_h=sqrt(tau_h/rohw)

!...total friction velocity in water
 	ustar_t=sqrt(ustar_w*ustar_w+ustar_h*ustar_h)
!
!...rain temperature = wet bulb temperature
 	tr=psyceq(p,ta,td)+t0

!...friction velocity related to kinetic energy of rain
 	ustar_r=rfric(rr)

!...friction velocity related to horizontal rain surface stress
!...coupling of total friction velocity with that by rain kinetic energy
 	ustar=c_ustar(ustar_t,ustar_h,ustar_r)

!...first guess of renewal time
 	tmax=t_star(q0,ql,qsol,ustar,tsurf,rr,sal)
 	call frain_ln(rr,tmax,sal,srs,srv)
 	dsal = srv + srs

!...compute delta t versus renewal time
 	call skin_ln(qsol,q0,tsurf,tmax,sal,tc,ts,gc)

!...calculate delta t versus rainrate
 	call rain_ln(tr,tsurf,rr,tmax,sal,trs,trv)

!...total temperature difference across the cool skin
	dt_skin=tc+ts+trs+trv

! second call to bulk-algorithm with parameterized skin temperature if
! skin_flag is false

	if(.not.skin_flag)then
		tsurf=tsurf+dt_skin
		esat=goff_gratch(tsurf)
		esat=esat*0.98
		qs=0.622*esat/(p-0.378*esat)
		call ren_bulk
	endif

! bubble-mediated gas-transfer velocity
	bc = bubble_k(ws, sal, tsurf)

 	end subroutine moblam8
!
!
!
       subroutine ren_bulk
!=========================
!	to evaluate surface fluxes, surface roughness and stability of
!	the atmospheric surface layer
!...version 11-jun-97, Peter Schluessel
!
!...Input: 
!	u, wind speed (m/s)
!	t, air temperature (K)
!	q, specific humidity (kg/kg)
!	ts, surface skin temperature (K)
!	qs, surface specific humidity (kg/kg)
!	rns, solar net flux (W/m^2)
!	rnl, longwave net flux (W/m^2)
!	rain, rain rate (mm/h)
!	sal, salinity (0/00)
!	zu, anemometer height (m)
!	zt, thermometer height (m)
!	zq, psychrometer height (m)
!	p, pressure (hPa)
!	skin_flag, if surface temperature is skin temperature
!...Output:
!	ustar_a: friction velocity of air (m/s)
!	tstar, MO scaling parameter for temperature (K)
!	qstar, MO scaling parameter for specific humidity (kg/kg)
!	zo, roughness length (m)
!	zot, roughness length for temperature (m)
!	zoq, roughness length for humidity (m)
!	zl, Obukhov stability parameter
!	tau, wind stress (N/m)
!	hf, sensible heat flux (W/m^2)
!	ef, latent heat flux (W/m^2)
!	tren, renewal time in air (s)
!
	real, parameter :: kar=0.40, t0=273.15, grav=9.806
	real, parameter :: rgas=287.05, zi=600.
    real, parameter :: beta=1.2! Fairall's low wind turbulence
	real kappa, epsil
	real obukhov, l, tau, hf, ef
	real st, stanton, da, dalton, tren, t_ren_air
	real lambda_t,lambda_v, sigma2_t, sigma2_v, a_t, a_v
	real eta_t, eta_v
	logical s_flag,skin_flag
	common/skin_in/u,t,q,ts,qs,rns,rnl,rain,sal,zu,zt,zq,p,skin_flag
	common/bulk_out/ustar_a,tstar,qstar,zo,zot,zoq,zl, tau, hf, ef,tren
	common/rparm/lambda_t,lambda_v,sigma2_t,sigma2_v,a_t,a_v,eta_t,eta_v,s_flag
!
!...set renewal parameters
	if(.not. s_flag) call ren_parm

	ta=t-t0
!--------------------------- air "constants" ---------------------------
!...latent heat of vaporization at skin temperature in J/kg
    hl=(2.5004-0.002363*(ts-t0))*1e+6

!...specific heat of dry air in J kg^-1 K^-1
	cpa=1005.

!...moist air density (kg/m^3)
	rhoa=p*100./(rgas*t*(1.+.608*q))

!...kinematic viscosity of dry air  
	visa = vis_air(t)

!...thermal diffusivity, m^2/s
	kappa=visa*1.4079

!...water vapour diffusivity in air	
	epsil=0.226e-4*(t/t0)**1.81*1013.25/p
!
!...get zo due to rain
 	zo_r=rain_z0(rain)

!...surface current
	us=0.

!...gustiness factor initial guess (Fairall et al., JGR, 101, 3747-3764, 1996)
    wg=0.5
    du=u-us
	du_wg=(du**2.+wg**2.)**.5

!...roughness initial guess
        zo=.0005

!...potential temperature diff
        dt=t-ts+.0098*zt
        dq=q-qs
        ustar_a=.04*du_wg

!...initial guesses
        tstar=.05*dt
        qstar=.05*dq


index_loop:       do index=1,20

!...zo at neutral conditions after Smith 1988
        zo=0.011*ustar_a*ustar_a/grav + 0.11*visa/ustar_a

!...modify if rain-induced roughness is higher
 	if(zo < zo_r)zo=zo_r

       zl=obukhov(t, q, ustar_a, tstar, qstar, zu)
       	if(abs(zl) > 1.e-10)then
 
! get stability functions
			l=zu/zl
			call stab_fun(dt, zo, zu, zt, l, psi_m, psi_h)
			psi_e=psi_h
       else
			psi_m=0.
			psi_h=0.
			psi_q=0.
       endif

!...gustiness effect (Fairall et al., JGR, 101, 3747-3764, 1996)
        ustar_a=du_wg*kar/((log(zu/zo)-psi_m))

! calculate renewal time in air
		tren=t_ren_air(t,tstar,ustar_a,lambda_t, a_t)
	
! caculate interfacial Stanton and Dalton numbers
		st=stanton(ustar_a, tren, kappa, eta_t)
		da=dalton(ustar_a, tren, epsil, eta_v)
		zot=zo*exp(kar*(5.0-1./(1.0*st)))
		zoq=zo*exp(kar*(5.0-1./(1.0*da)))
		s = (log(zt/zot)-psi_h)/kar
		d = (log(zq/zoq)-psi_e)/kar
 
!... modify fluxes
       tstar=dt/s
       qstar=dq/d
       tvsr=tstar*(1.+0.608*q)+(0.608*t*qstar)
       bf=-grav/ta*ustar_a*tvsr

!...calculate gustiness factor
        if(bf.gt.0) then
			wg=beta*(bf*zi)**0.333
        else
			wg=0.
        endif
!...include gustiness in wind spd.
        du_wg=(du**2.+wg**2.)**.5

       enddo index_loop

!...calculate fluxes
            tau=rhoa*ustar_a*ustar_a
            hf=cpa*rhoa*ustar_a*tstar
            ef=hl*rhoa*ustar_a*qstar

       end subroutine ren_bulk 
!
!
!
      real function obukhov(ta, q, ustar, tstar, qstar, z) result (zl)
!=====================================================================
! calculates Obukhov stability parameter z/l
! all quantities in SI units
! ta = air temperature
! q = specific humidity
! z = height
! ustar = friction velocity
! tstar = frictional temperature
! qstar = frictional moisture
!       see Liu et al. (1979)

	implicit none
	real ta, q, ustar, tstar, qstar, z
	real, parameter :: t0 = 273.15, grav = 9.806, kar = 0.4
	real tv, tvstar, ob


      tv=ta*(1.+0.608*q)
      tvstar=tstar*(1.+0.608*q)+0.608*ta*qstar

      if(tvstar /= 0. )then
		ob=tv*ustar*ustar/(grav*kar*tvstar)
		zl=z/ob
      else
		zl = 0.
      endif

      end function obukhov
!
!
!
	subroutine stab_fun(dt, z0, zm, zh, l, psi_m, psi_h)
!=======================================================
! dt = t_skin - t_air
! z0 = roughness length
! zm = anemometer height
! zh = psychrometer height
! l  = Monin-Obukhov length
! psi_m = stability function for momentum
! psi_h = stability function for heat and moisture
!
	implicit none
	real dt, z0, zm, zh, l, psi_m, psi_h
	real, parameter :: a=0.7, b=0.75, c=5., d=0.35
	real, parameter :: ah=1.0, bh=0.667
	real zw, zwm, zwh, zw0, zw1, zeta, zeta0, lambda, lambda0

	if( dt < 0.) then

! unstable condition, Benoit (1977)

	zwm=1.-15.*zm/l
	zw0=1.-15.*z0/l


	zeta = zwm**0.25
	zeta0 = zw0**0.25

	psi_m =  -log((zeta0*zeta0+1)*(zeta0+1)*(zeta0+1)/ &
	((zeta*zeta+1)*(zeta+1)*(zeta+1)))                &
	 -2.*(atan(zeta)-atan(zeta0))

	zwh=1.-9.*zh/l
	zw0=1.-9.*z0/l
	lambda = zwh**0.5
	lambda0 = zw0**0.5

	psi_h = -2.*log((lambda0+1)/(lambda+1.))

	elseif (dt > 0.) then

! stable condition, Beljaars and Holtslag (1991)

	zw = zm/l
	zw1 = b*c/d

	zw1 = -(a*zw+b*(zw-c/d)*exp(-d*zw)+zw1)
	zw=zm/z0
	psi_m=log(zw/(zw-zw1))
	zw = zh/l
	zw1=(1.+2.*ah*zw/3.)**(2/3)
	zw1 = -(zw+bh*(zw-c/d)*exp(-d*zw)+zw1-1.)
	zw=zh/z0
	psi_h=log(zw/(zw-zw1))

	else

! neutral condition
	
	psi_m = 0.
	psi_h = 0.
	
	endif
	
	end subroutine stab_fun
!
!
!
 	real function rain_z0(rr)
!============================
!...calculation of roughness lenght for rain-induced waves
!...input: rainrate in mm/h
!...output: roughness lenght in m
!
 	real, parameter:: ex=-0.21d0, rc=0.30d0, ups=1.15045d0
 	real, parameter:: a0=0.0129d0, a1=1.60686d0, a2=-1.76407d0
 	real, parameter:: dvier=3.d0/4.d0
	real, parameter :: pi=3.14159265358979
 	external gamma_ln

	if(rr < 0.01)then
	rain_z0 = 0.
	return
	endif

	svier=7.d0/4.d0
 	rate=rr
 	gam=4.1d0*rate**ex
 	fac1=2.d0*gam/exp(-2.d0*gam*rc)
 	arg1=2.d0*gam*rc
 	arg2=(ups+2.d0*gam)*rc
 	fac2=a1*(2.d0*gam)**(-svier)
 	fac3=a2*(ups+2.d0*gam)**(-svier)
 	hwred=a0+fac1*(fac2*gamma_inc(svier,arg1)+fac3*gamma_inc(svier,arg2))

!...integration over all drop radii
 	z0red=0.058*hwred**1.19d0
 	gsvier=exp(gamma_ln(svier))

!...integration over drops with radii greater than r_c
 	z0redv=a0+gsvier*(a1*(2.d0*gam)**(-dvier)+                             &
       	2.d0*a2*gam*(ups+2.d0*gam)**(-svier))
!
!	rain_z0=z0redv*1.d-3
 	rain_z0=z0red*1.d-3

 	end function rain_z0
!
!
!
       real function gamma_ln(xx)
!================================
!...supports function gamma_inc
       real xx
       integer j
       real ser,tmp,x,y
       real,parameter::stp=2.5066282746310005d0
       real,parameter,dimension(6):: cof=(/76.18009172947146d0,-86.50532032941677d0,&
       24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, -.5395239384953d-5/)

       x=xx
       y=x
       tmp=x+5.5d0
       tmp=(x+0.5d0)*log(tmp)-tmp
       ser=1.000000000190015d0
       do j=1,6
       y=y+1.d0
       ser=ser+cof(j)/y
       enddo
       gamma_ln=tmp+log(stp*ser/x)

       end function gamma_ln
!
!
!
 	real function gamma_inc(a,x)
!===============================
!...incomplete gamma function of real order

 	external gamma_ln,gamma_q

 	zw=gamma_ln(a)
 	gamma_inc=gamma_q(a,x)*exp(zw)

 	end function gamma_inc
!
!
!
       real function gamma_q(a,x)
!================================
!...supports function gamma_inc

       real a,x
       external hmg,gseries
       real gammcf,gamser,gln

       if(x < a+1.)then
         call gseries(gamser,a,x,gln)
         gamma_q=1.-gamser
       else
         call hmg(gammcf,a,x,gln)
         gamma_q=gammcf
       endif

       end function gamma_q
!
!
!
       subroutine hmg(gammcf,a,x,gln)
!====================================
!...supports function gamma_q

       real a,gammcf,gln,x,eps,fpmin
       external gamma_ln
       integer, parameter:: itmax=1000
       integer i
       real an,b,c,d,del,h,gamma_ln

		eps=1.d-15
		fpmin=tiny(x)
		gln=gamma_ln(a)
		b=x+1.-a
		c=1./fpmin
		d=1./b
		h=d

       do i=1,itmax
		an=-i*(i-a)
		b=b+2.
		d=an*d+b
       if(abs(d) < fpmin)d=fpmin
		c=b+an/c
       if(abs(c) < fpmin)c=fpmin
		d=1./d
		del=d*c
		h=h*del
       if(abs(del-1.).lt.eps)exit
       enddo
	
		gammcf=exp(-x+a*log(x)-gln)*h

       end subroutine hmg
!
!
!
       subroutine gseries(gamser,a,x,gln)
!========================================
!...supports function gamma_q

       real a,gamser,gln,x,eps
       integer, parameter:: itmax=1000
       external gamma_ln
       integer n
       real ap,del,sum,gamma_ln

       eps=1.d-15

       gln=gamma_ln(a)

       if(x <= 0.)then
       if(x < 0.)pause 'x < 0 in gseries'
       gamser=0.
       return
       endif

       ap=a
       sum=1./a
       del=sum

       do n=1,itmax
		ap=ap+1.
		del=del*x/ap
		sum=sum+del
		if(abs(del) < abs(sum)*eps)exit
       enddo

       gamser=sum*exp(-x+a*log(x)-gln)

       end subroutine gseries
!
!
!
	function t_ren_air(t,tstar,ustar, lambda_0, a) result(tmax)
!==============================================================
!...calculates renewal time in air
!
	real, parameter :: rfcr = -1.5e-4, zw=27./128.
	real lambda_0, a
	real, parameter :: grav = 9.806, pi = 3.14159265358979
	real rf0, alfa, nue, zw1

    nue = vis_air(t)
	alfa=1./t
	rf0=alfa*grav*tstar*nue/ustar**3.
	zw1=max(1.e-40,1.+rf0/rfcr)

	tmax=zw*a*nue*pi*pi*lambda_0*lambda_0/ &
		(sqrt(zw1)*(ustar*ustar))
	tmax=min(tmax,10000.)

	end function t_ren_air
!
!
!
	function stanton(ustar, tren, kappa, eta) result(sta)
!========================================================
!...calculate Stanton number

	real ustar, tren, sta
	real kappa, eta

	sta=eta*sqrt(kappa/(ustar*ustar*tren))

	end function stanton
!
!
!
	function dalton(ustar, tren, epsil, eta) result(dal)
!=======================================================
!...calculate Dalton number

	real ustar, tren
	real epsil, eta

	dal=eta*sqrt(epsil/(ustar*ustar*tren))

	end function dalton
!
!
!
 	subroutine  skin_ln(qsol,q0,temp,tmax,sal,tc,ts,gas_coef)
!================================================================
!...calculate cool and warm skin
!...A. V. Soloviev (Matlab Version)
!...P. Schluessel (Fortran)
!
 	real, parameter :: t0=273.15d0, pi=3.14159265358979
 	real,parameter,dimension(9)::a3 = (/0.237,0.360,0.179,0.087,0.080,&
 		0.0246,0.025,0.007,0.0004/)
 	real,parameter,dimension(9)::alfa3 = (/2.8735e-2,4.40529e-1,3.1746e+1,&
 		182.48,1201.92,7936.5,3194.89,12788.,69444./)
 	real k, mu
 	external erf
!
!...density of sea water in kg/m^3
 	roh=densit(temp,sal)

!...specific heat of sea water in J/(kg K)
 	cp = spec_h(temp,sal)

!...molecular gas diffusion coefficient for carbon dioxide in m^2/s
	mu = diff_co2(temp)

!	print *, "gas diffusivity in water: ", mu
!...normalized heat fluxes
 	qr=qsol/cp/roh
 	q=q0/cp/roh

!...thermal diffusion coefficient
         k=(0.1344e-2+0.42e-5*(temp-t0))*1.e-4

!...cool skin without solar heating but with modified Rf0
 	sigma2=6.23
 	tc=q*4./3.*(tmax*exp(-sigma2/16.)/k)**0.5*pi**(-0.5)

!...gas-transfer coefficient
	gas_coef=2.*sqrt(pi)*exp(3.*sigma2/16.)*sqrt(mu/tmax)

!...calculation of solar haeting with log-normal weighting
 	tsmean=0.

 	do jj=1,20
 	v=-3.+(jj+1)/20.*6.
 	tr=tmax*exp(sqrt(sigma2)*v-sigma2/4.)
 	tsm=0.
 	 do ii=1,9
 	  del=alfa3(ii)*(k*tr)**0.5
 	  delq=del*del
 	  if(del > 3)then
 		tsm=tsm+qr*a3(ii)/alfa3(ii)/k*((1./sqrt(pi)/del*                 &
       	(1.-1./(2.*delq)+3./(2.*delq)**2.-15./(2.*delq)**3.)             &
       	-1.)/(delq)-1.+2./sqrt(pi)/del+4./3./sqrt(pi)*del)
 	  else
 		if(del > 0.01)then
 			tsm=tsm+qr*a3(ii)/alfa3(ii)/k*((exp(delq)-1.)/delq               &
       		-1.+(2./sqrt(pi)-exp(delq)*erf(del)/del)/del+4./3./sqrt(pi)      &
       		*del)
 		else
 			tsm=tsm+qr*a3(ii)/alfa3(ii)/k*(delq)/2.
 		endif
 	  endif
 	 enddo
 	  tsmean=tsmean+0.3/sqrt(pi)*exp(-v*v)*tsm
 	enddo

 	ts=tsmean

 	end subroutine skin_ln
!
!
!
 	real function t_star(q0,ql,qsol,ustar,t,rr,sal)
!=====================================================
!...compute renewal time for given heat flux, friction velocity,
!...and temperature
!...input: q0 = heat flux in w/m^2, upward negative
!...       ql = latent heat flux in w/m^2, upward negative
!...     qsol = solar heat flux in w/m^2, downward positive
!...       ustar = friction velocity in upper ocean in m/s
!...       t = surface temperature in k
!...       rr = rainrate in mm/h
!...       sal = salinity in parts per thousand
!...output: renewal time in s
!...Peter Schluessel
!...20-feb-96
!
 	implicit none
 	real, parameter:: bettas=0.00075, t0=273.15
 	real, parameter:: pi=3.14159265358979, g=9.806, lam0=13.3
 	real, parameter:: kecr=0.18d0, rfcr=-1.5d-4
 	real, parameter:: sigma2=6.23
 	real q0,ustar,t,k,ql,qeff,qsol,xlat,ql0,ra,zmax, p
 	real ke,alphat,xnue,cp,roh,qq0,rf0
 	real sal,densit,viscos,vis,d,tstarq,tstarr,rr,rain
 	real txrain,spec_h
 	external txrain, densit, viscos, spec_h, cdepth
!...freezing point
!	t0=273.15
!...coefficient of salinity expansion
!	bettas=0.00075
!...pressure in dbar
	p=10.
!...gravity in m/s^2
!	g=9.806
!	lam0=13.3
!...kritical Keulegan number
!	kecr=0.18
!...kritical surface Richardson number
!	rfcr=-1.5e-4
!-----------------------------------------------------------------------

!...rainrate in m/s
 	rain=rr/3.6d6

!...molecular heat transfer coefficient in m^2/s
         k=(0.1344e-2+0.42e-5*(t-t0))*1.e-4

!...latent heat of evaporation
 	xlat=2.5004e6-2.3634e3*(t-t0)

!...density of sea water in kg/m^3
 	roh=densit(t,sal)

!...dynamic viscosity in kg/(m s)
 	vis=viscos(t,sal,p)

!...kinematic viscosity in m^2/s
 	xnue=vis/roh

!...coefficient of thermal expansion of sea water in K^-1
 	alphat=(0.0977*(t-t0)+0.52)*1.e-4

!...specific heat of sea water in j/(kg k)
 	cp=spec_h(t,sal)

!...compensation depth and effective heat flux
 	if(qsol > q0 .and. q0+qsol > 2.)then
  		call cdepth(q0,qsol,d,zmax,qeff)

!...Rayleigh number
 		ra=alphat*g*(-qeff/cp/roh)*zmax**4./xnue/k**2.

 		if(ra < 1300)then
 			qq0=0.
 		else

!...normalized heat flux
 			qq0=-q0/(cp*roh)
 		endif
 	else
 		qq0=-q0/(cp*roh)
 	endif

!...normalized latent heat flux
 	ql0=-ql/(xlat*roh)

!...surface Richardson number
 	rf0=(-alphat*qq0-bettas*sal*(ql0-rain))*g*xnue/ustar**4.
 	rf0=min(rf0,0.)

!...Keulegan number
 	ke=ustar**3./(g*xnue)

!...renewal time in s
 	tstarq=9./16.*exp(sigma2/8.)*xnue*pi*lam0*lam0*          &
       	(1+rf0/rfcr)**(-0.5)*(1+ke/kecr)/(ustar*ustar)

!...get renewal time and combine both surface renewals
 	if(rr > 1.d-5)then
 		tstarr=txrain(rr)
 		t_star=1.d0/(1.d0/tstarq+1.d0/tstarr)
 	else
 		t_star=tstarq
 	endif

 	end function t_star
!
!
!	
 	subroutine rain_ln(tr,ts,rr,tmax,sal,trs,trv)
!====================================================
!...surface and volume heat flux due to rain
!...P. Schluessel
!...20-feb-96
!
!
!...input: tr = rain temperature in k
!          ts = surface temperature in k
!          rr = rainrate in mm/h
!          tmax = renewal time in s
!...output: trs = skin temperature deviation due to surface cooling in k
!           trv = skin temperature deviation due to volume cooling in k
!
 	save
 	real, parameter:: rc=0.4d0, ex=-0.21d0, t0=273.15d0, a=1.d2
 	real, parameter:: sigma2=6.23d0, pi=3.14159265358979

!...zeta in mm^-1
        real, parameter,dimension(14):: zeta =(/-0.67227554e-04,-0.12280522e-02,&
        -0.77123456e-03,-0.57451221e-03,-0.22742890e-02,-0.64955258e-04,&
        -0.87949847e-03,-0.40283216e-06,-0.27432975e-05,-0.15406247e-02,&
       -0.43799929e-03,-0.26090549e-02,-0.44233906e-04,-0.33905604e-04/)
        real, parameter,dimension(14):: psi =(/0.59670148e+00,0.18085820e+01,&
        0.33960565e+00,0.96367711e+00,0.29968117e+00,0.27661104e+01,&
        0.21749846e+00,0.65308387e+01,0.26095015e+01,0.55630270e+00, &
       0.65080732e+00,0.60107111e+00,0.22811756e+01,0.27673898e+01/)
 	real k
 
!
!...rainrate in m/s
 	rain=rr*1.d-3/3.6d3
 	trs=0.d0
 	trv=0.d0

!...gam in mm^-1
 	gam=4.1d0*rr**ex
 	gamrc=gam*rc

!...kappa in m^2/s
 	k=(0.1344d-2+0.42d-5*(ts-t0))*1.d-4

!...density of sea water in kg/m^3
 	roh=densit(ts,sal)

!...density of fresh water in kg/m^3
 	rohf=1.d3

!...specific heat of sea water in j/(kg k)
 	cp=spec_h(ts,sal)

!...specific heat of fresh water in j/(kg k)
 	cpf=spec_h(ts,0.d0)

!...fresh versus sea water ratio
 	fsratio=cpf*rohf/(cp*roh)

!...surface heat flux due to rain
 	trs=rain*(tr-ts)*fsratio*4.d0/3.d0*sqrt(tmax/(k*pi))*                  &
       	(1.d0+2.d0*gamrc+2.d0*gamrc*gamrc+4.d0*gamrc*gamrc*gamrc/3.d0)*  &
       	exp(-2.d0*gamrc)*exp(-sigma2/8.d0)

!...gam in m^-1
 	gam=gam*1.d3
 
 	trmean=0.d0
         do jj=1,20
         v=-3.+(jj+1)/20.*6.
         tx=tmax*exp(sqrt(sigma2)*v-sigma2/4.)
         tsm=0.
         do ii=1,14
 	del=psi(ii)*gam*sqrt(k*tx)/a
 	delq=del*del
         if(del > 3.)then
         tsm=tsm+(1.d3*zeta(ii)*a/(gam*psi(ii)*k))*((1./sqrt(pi)/del*    &
         (1.-1./(2.*delq)+3./(2.*delq)**2.-15./(2.*delq)**3.)            &
         -1.)/(delq)-1.+2./sqrt(pi)/del+4./3./sqrt(pi)*del)
         else
         if(del > 0.01)then
         tsm=tsm+(1.d3*zeta(ii)*a/(gam*psi(ii)*k))*((exp(delq)-1.)/delq  &
         -1.+(2./sqrt(pi)-exp(delq)*erf(del)/del)/del+4./3./sqrt(pi)     &
         *del)
         else
         tsm=tsm+(1.d3*zeta(ii)*a/(gam*psi(ii)*k))*(delq)/2.
         endif
         endif
         enddo
 
 	trmean=trmean+0.4d0*exp(-v*v)*tsm/sqrt(pi)
 	enddo
 	trv=-trmean*rain*(tr-ts)*fsratio

 	end subroutine rain_ln
!
!
!
 	subroutine frain_ln(rr,tmax,sal,srs,srv)
!=====================================================
!...surface and volume fresh-water flux due to rain
!...P. Schluessel
!...20-feb-96
!
!
!...input: rr = rainrate in mm/h
!          tmax = renewal time in s
!          sal = salinity in parts per thousand
!...output: srs = salinity deviation due to surface flux
!           srv = salinity deviation due to volume flux
!
 	real, parameter:: rc=0.4d0, ex=-0.21d0, t0=273.15d0, a=1.d2
 	real, parameter:: pi=3.14159265358979, sigma2=6.23d0

!...zeta in mm^-1
        real, parameter, dimension(14):: zeta = (/-0.67227554e-04,-0.12280522e-02,&
        -0.77123456e-03,-0.57451221e-03,-0.22742890e-02,-0.64955258e-04,&
        -0.87949847e-03,-0.40283216e-06,-0.27432975e-05,-0.15406247e-02,&
       -0.43799929e-03,-0.26090549e-02,-0.44233906e-04,-0.33905604e-04/)
        real, parameter, dimension(14):: psi = (/0.59670148e+00,0.18085820e+01,&
        0.33960565e+00,0.96367711e+00,0.29968117e+00,0.27661104e+01,0.21749846e+00,&
       0.65308387e+01,0.26095015e+01,0.55630270e+00,0.65080732e+00,0.60107111e+00,&
       0.22811756e+01,0.27673898e+01/)
!
!
!...rainrate in m/s
 	rain=rr*1.d-3/3.6d3
 	srs=0.d0
 	srv=0.d0

!...gam in mm^-1
 	gam=4.1d0*rr**ex

!...gam in m^-1
 	gamrc=gam*rc

!...sigma in m^2/s
 	sig=1.5d-9

!...surface salinity flux due to rain
 	srs=-rain*sal*4.d0/3.d0*sqrt(tmax/(sig*pi))*                           &
       	(1.d0+2.d0*gamrc+2.d0*gamrc*gamrc+4.d0*gamrc*gamrc*gamrc/3.d0)*  &
       	exp(-2.d0*gamrc)*exp(-sigma2/8.d0)

!...gam in m^-1
 	gam=gam*1.d3
 
 	srmean=0.d0
         do jj=1,20
         v=-3.+(jj+1)/20.*6.
         tx=tmax*exp(sqrt(sigma2)*v-sigma2/4.)
         ssm=0.
         do ii=1,14
 	del=psi(ii)*gam*sqrt(sig*tx)/a
 	delq=del*del
         if(del > 3)then
         ssm=ssm+(1.d3*zeta(ii)*a/(gam*psi(ii)*sig))*((1./sqrt(pi)/del*  &
         (1.-1./(2.*delq)+3./(2.*delq)**2.-15./(2.*delq)**3.)            &
         -1.)/(delq)-1.+2./sqrt(pi)/del+4./3./sqrt(pi)*del)
         else
         if(del > 0.01)then
         ssm=ssm+(1.d3*zeta(ii)*a/(gam*psi(ii)*sig))*((exp(delq)-1.)/delq&
         -1.+(2./sqrt(pi)-exp(delq)*erf(del)/del)/del+4./3./sqrt(pi)     &
         *del)
         else
         ssm=ssm+(1.d3*zeta(ii)*a/(gam*psi(ii)*sig))*(delq)/2.
         endif
         endif
         enddo
 
 	srmean=srmean+0.4d0*exp(-v*v)*ssm/sqrt(pi)
 	enddo
 	srv=sal*srmean*rain

 	end subroutine frain_ln
!
!
!
      function goff_gratch(t) result(eswat)
!==========================================
!...saturation vapour pressure over water
!...Goff-Gratch formula
!... t = temperature (K)
!... eswat = saturation vapour press. over water (hPa)
!... Smithsonian Meteorological Tables,  p. 350

      real, parameter:: c1=7.90298, c2=5.02808, c3=1.3816e-7, c4=11.344
      real, parameter:: c5=8.1328e-3, c6=3.49149
!
      rmixv = 373.16/t
      es = -c1* (rmixv-1.) + c2*log10(rmixv) -&
          c3* (10.** (c4* (1.-1./rmixv))-1.) +&
          c5* (10.** (-c6* (rmixv-1.))-1.)
      eswat = 1013.246*10.**es
      return
      end function goff_gratch
!
!
!
 	function txrain(rate)
!============================
!...calculation of renewal time due to rain
!...P. Schluessel
!...16-mar-98

 	save
 	real, parameter:: pi=3.14159265358979, fac=4.1d0, ex=-0.21d0
 	real, parameter:: grav=9.806d0, alpha=9430.d0
 	real, parameter:: rc=0.40d0, ups=1.15045, a=1.0117d0, b=1.05201d0
 	real, parameter:: phi=1.278d0
 	real, parameter:: n0=16.e-6
 	integer, parameter:: n3=3, n4=4
 	external gamma_inc,gammain
	fuha = 5.d0/2.d0
 	sgrav=sqrt(grav*1.e3)
 	gam=fac*rate**ex
 	arg=2.d0*gam*rc
 	fac1=sgrav/(pi*n0*phi*phi*alpha*alpha)
 	sum1=a*a*(2.d0*gam)**(-fuha)*gamma_inc(fuha,arg)
 	arg=(2.d0*gam+ups)*rc
 	sum2=-2*a*b*(2.d0*gam+ups)**(-fuha)*gamma_inc(fuha,arg)
 	arg=(2.d0*gam+2.d0*ups)*rc
 	sum3=b*b*(2.d0*gam+2.d0*ups)**(-fuha)*gamma_inc(fuha,arg)

 	txrain=fac1/(sum1+sum2+sum3)

 	end function txrain
!
!
!
       function erf(x)
!=====================
!...calculates error function
!
       real erf,x
       real gamma_p
       external gamma_p
       if(x < 0.)then
       erf=-gamma_p(.5d0,x**2.d0)
       else
       erf=gamma_p(.5d0,x**2.d0)
       endif

       end function erf
!
!
!
       function gamma_p(a,x)
!===========================
       real a,gamma_p,x
       external gseries,hmg
       real gammcf,gamser,gln
       if(x.lt.0..or.a.le.0.)pause 'bad arguments in gamma_p'
       if(x.lt.a+1.)then
       call gseries(gamser,a,x,gln)
       gamma_p=gamser
       else
       call hmg(gammcf,a,x,gln)
       gamma_p=1.-gammcf
       endif

       end function gamma_p
!
!
!
       function psyceq(p,t,td) result(ewt)
!=========================================
!...calculation of wet bulb temperature from temperature (deg c),
!...dewpoint (deg c), and pressure (hpa) using the psychrometric
!...equation
!
         tmin = td
 	diff=9999999.d0
 	do while(diff >= 0.1)
         ewt = tmin + 0.01
         e = 6.11d0 * 10.d0**((7.5d0 * td)/(237.7d0 + td))
         esw = 6.11d0 * 10.d0**((7.5d0 * ewt)/(237.7d0 + ewt))
         edif = esw - e
         peq = 0.000660d0 * (1.d0 + 0.00155d0 * ewt) * p * (t - ewt)
         diff = abs(peq - edif)
         tmin = ewt
 	enddo

         end function psyceq
!
!
!
 	function densit(tk,s)
!============================
!...temperature in k
!...salinity in per mille (parts per thousand)
!...output: density of sea water in kg/m^3 at atmospheric pressure
!...from Dietrich, Kalle, Krauss, Siedler

	real, parameter:: t0=273.15
 	real ai(4),aij(2,0:3),b(0:3)
 	t=tk-t0
	b=0.
 	b(0)=-0.0934458632
 	b(1)=0.814876577
 	b(2)=-4.82496140e-4
 	b(3)=6.76786136e-6
 	sigma0=0.
 	do i=0,3
    	sigma0=sigma0+b(i)*s**float(i)
    	enddo
!
	ai=0.
	aij=0.
 	ai(1)=4.53168426
 	ai(2)=-0.54593911
 	ai(3)=-1.98248399e-3
 	ai(4)=-1.43803061e-7
 	aij(1,0)=1.
 	aij(1,1)=-4.7867e-3
 	aij(1,2)=9.8185e-5
 	aij(1,3)=-1.0843e-6
 	aij(2,1)=1.8030e-5
 	aij(2,2)=-8.164e-7
 	aij(2,3)=1.667e-8
 	x=0.
 	do i=1,4
   	x=x+ai(i)*t**float(i)
   	enddo
 	x=x/(t+67.26)
 	y=0.
 	do i=1,2
 	do j=0,3
   	y=y+aij(i,j)*sigma0**float(i)*t**float(j)
   	enddo
   	enddo
 	sigmat=x+y
 	densit=(1.+sigmat*0.001)*1000.

 	end function densit
!
!
!
 	real function spec_h(t,sal) result(cp0t)
!===============================================
!...specific heat of fresh and sea water at surface pressure
!...Millero et al. (1973)
!...t in K
!...sal in parts per thousand
!
 	save
 	real, parameter:: t0=273.15
 	tt=t-t0
 	tt2=tt*tt
 	tt3=tt2*tt
 	tt4=tt3*tt
 	cp0t=4217.4d0-3.720283*tt+0.1412855d0*tt2-2.654387d-3*tt3+             &
       	2.093236d-5*tt4
 	if(sal > 0.d0)then
 	cp0t=cp0t+sal*(-7.644d0+0.107276*tt-1.3839d-3*tt2)+                  &
       	sal**1.5d0*(0.17709d0-4.0772d-3*tt+5.3539d-5*tt2)
 	endif

 	end function spec_h
!
!
!
 	real function viscos(tk,s,p)
!===================================
!...temperature in K
!...salinity in per mille
!...pressure in dbar
!...output: dynamic viscosity of sea water in kg/(m s)
!
 	t=tk-273.15
 	t2=t*t
 	t3=t2*t
 	p2=p*p
 	v=1.79e-2-6.1299e-4*t+1.4467e-5*t2-1.6826e-7*t3                        &
       	-1.8266e-7*p+9.8972e-12*p2+2.4727e-5*s+                          &
         s*(4.8429e-7*t-4.7172e-8*t2+7.5986e-10*t3)+                     &
         p*(1.3817e-8*t-2.6363e-10*t2)-                                  &
       	p2*(6.3255e-13*t-1.2116e-14*t2)
 	viscos=v*0.1

 	end function viscos
!
!
!
  	subroutine cdepth(q0,qsol,d,zmax,qeff)
!=============================================
!...calculates compensation depth and effective heat flux
!...by Newton's iteration
!...input parameters: q0 - surface heat flux due to longwave radiation,
!                         sensible and latent heat flux (in W/m^2)
!                    qsol - solar radiation flux penetrating ocean
!                           surface  (in W/m^2)
!...output parameters: d - compensation depth (in m)
!                      qeff - effective heat flux  (in W/m^2)
!                      zmax - depth of maximum of Rayleigh number (m)
!
!
 	real, parameter, dimension (9):: a3 = (/0.237,0.360,0.179,0.087,0.080,&
 	0.0246,0.025,0.007,0.0004/)
 	real, parameter, dimension(9):: alfa3 = (/2.8735e-2,4.40529e-1,3.1746e+1,&
 	182.48,1201.92,7936.5,3194.89,12788.,69444./)


 	q=q0
 	qr=qsol
!...compensation depth
 	d=0.

 	do 
 	sumf=0.
 	sumfg=0.
 	f=qr*(1.-sum(a3*exp(-alfa3*d)))+q
 	fg=qr*sum(a3*alfa3*exp(-alfa3*d))
 	r=f/fg
 	d=d-r
 	if(abs(r/d) < 0.001)exit
 	enddo
!
 	z=d
 	do
 	sum1=sum(a3*exp(-alfa3*z)-exp(-alfa3*d))
 	sum2=sum(a3*alfa3*exp(-alfa3*z))
 	sum3=sum(alfa3*a3*alfa3*exp(-alfa3*z))
 	f=4.*z**3.*sum1-z**4.*sum2
 	fg=12.*z**2.*sum1-4.*z**3.*sum2-4.*z**3.*sum2+z**4.*sum3
 	r=f/fg
 	z=z-r
 	if(abs(r/z) < 0.01)exit
 	enddo
 	zmax=z
!
 	qeff=qr*sum(a3*(exp(-alfa3*d)-exp(-alfa3*zmax)))

 	end subroutine cdepth
!
!
!
 	function rfric(rr)
!===========================
!...calculate friction velocity related to the kinetic energy of the
!...rain
!...input: rain rate in mm/h
!...output: friction velocity in m/s
!
!...P. Schluessel, 26-mar-98
!
 	real, parameter:: pi=3.14159265358979
 	real, parameter:: rc=0.40d0, ex=-0.21d0,xn0=16.d3
 	real, parameter:: b1=1.0117d0, b2=1.05201d0, ups=1.15045d0
 	real, parameter:: alpha=9.43

!...alpha in m/s, ups in 1/mm, rc in mm
!...gam in mm^-1
         gam=4.1d0*rr**ex
         gamrc=gam*rc
 	zgam=2.d0*gam

!...gam in m^-1
 	gam=gam*1.d3
!
 	fac1=2.d0*pi*xn0*alpha*alpha*alpha/3.d0
 	sum1=b1*b1*b1*gammain(4,2.d0*gamrc)/(16.d0*gam**4.d0)
 	sum2=-3.d0*b1*b1*b2*gammain(4,(zgam+ups)*rc)/                             &
       	(2.d0*gam+ups*1.d3)**4.d0
 	sum3=3.*b1*b2*b2*gammain(4,(zgam+2.d0*ups)*rc)/                              &
       	(2.d0*gam+2.d3*ups)**4.d0
       	sum4=-b2*b2*b2*gammain(4,(zgam+3.d0*ups)*rc)/            &
       	(2.d0*gam+2.d3*ups)**4.d0
 	xke=fac1*(sum1+sum2+sum3+sum4)
 	rfric=(xke*1.d3)**(1.d0/3.d0)

 	end function rfric
!
!
!
 	function gammain(np1,x)
!==============================
!...incomplete gamma-function of integer order
 	n=np1-1
 	xexp=exp(-x)
 	fakult=1.d0
 	xm=1.d0
 	sum=0.d0
 	do i=0,n
 	if(i > 0)then
 	fakult=fakult*float(i)
 	xm=xm*x
 	endif
 	sum=sum+xm/fakult
 	enddo
 	gammain=fakult*xexp*sum

 	end function gammain
!
!
!
 	function rstress(u10,r)
!==========================
!...calculates surface stress produced by rainfall
!...input: u10: wind speed at 10 m in m/s
!          r: rain-rate in mm/h
!...output: rstress in N/m^2
! P. Schluessel, 7-mar-96
!
 	real, parameter:: gamma=0.85d0, rho=1.d3
 	rstress=gamma*rho*u10*r*1.d-3/3600.d0

 	end function rstress
!
!
!
 	function c_ustar(ustar_t,ustar_h,ustar_r)                                      
!============================================
!...routine couples friction velocities
!...ustar_t: wind-induced friction velocity
!...ustar_h: friction velocity related to horizontal momentum of rain 
!...ustar_r: friction velocity related to kinetic energy of vertical
!...         rain velocity
!...c_ustar: coupled friction velocity
!
 	parameter (drittel = 1.d0/3.d0)                                        
!...enhanced wave breaking is only considered if wind-induced
!...friction velocity exceeds certain threshold: 0.5 cm/s assumed
!
 	thresh=0.005                                                           
!
 	if(ustar_t.gt.thresh)then                                              
 	c_ustar=(ustar_t**3.d0+ustar_h**3.d0+ustar_r**3.d0)**drittel
 	else                                                                   
 	c_ustar=(ustar_t**3.d0+ustar_h**3.d0)**drittel                         
 	endif                                                                  
!
 	return                                                                 
 	end                                                                    
!
!
!
	function dewpt(q,p) result(t)
!================================
! input q in kg/kg
!       p in hPa
! output dewpoint in K
!
	e=q*p/(0.622+0.378*q)
	z=log10(e/6.107)
	t=237.*z/(7.5-z)
	t=t+273.15
	return
	end function dewpt
!
!
!
	subroutine ren_parm
!======================

	common/rparm/lambda_t,lambda_v,sigma2_t,sigma2_v,a_t,a_v,eta_t,eta_v,s_flag
	real lambda_t,lambda_v, sigma2_t, sigma2_v, a_t, a_v
	real eta_t, eta_v
	logical s_flag
	real, parameter:: cr=122., pi=3.14159265358979
	real, parameter:: rfcr=-0.5524

	lambda_t=(-rfcr)**(-0.25)*0.144**(-0.75)
	lambda_v=(-rfcr)**(-0.25)*0.159**(-0.75)
	

	sigma2_t=-16.*log(lambda_t*0.75*sqrt(pi/cr))   
	sigma2_v=-16.*log(lambda_v*0.75*sqrt(pi/cr))

	a_t=8.*exp(sigma2_t/8)/(pi*3.)
	a_v=8.*exp(sigma2_v/8)/(pi*3.)

	eta_t=2.*exp(sigma2_t*3./16.)
	eta_v=2.*exp(sigma2_v*3./16.)

	s_flag=.true.

	end subroutine ren_parm
!
!
!
	real function bubble_k (wind, sal, temp)
!===========================================
! bubble mediated gas transfer velocity according to
! Woolf, D.K. (1997): Bubbles and their role in gas exchange, in: The sea surface
! and global change, ed. by P.S. liss and R.A. Duce, Cambridge University Press, Cambridge,
! 173-205
! input wind speed in m/s, salinity in psu, temperature in K

	real wind, sal, temp
	real, parameter:: v_factor = 6.25e-3, e = 1.4e4, f = 1.2
	real, parameter:: clean_n = 0.5, dirty_n = 0.66666667
	real, parameter:: r_gas = 8.314
	real solubility, v, diffusivity, whitecap, foam, solub, sol_co2, diff_co2, n_exp
!
! assuming clean or dirty bubbles
	n_exp = dirty_n
!
! whitecap coverage
	whitecap = foam(wind)
!
! bubble volume flux density
	v = v_factor * whitecap

! solubility in mol/(m^3 Pa)
	solub = sol_co2(sal, temp)

! Ostwald solubility  (dimensionless)
	solubility = solub * r_gas * temp

! diffusion coefficient
	diffusivity = diff_co2(temp)
! 
	bubble_k = v * (1. + (e * solubility * diffusivity ** n_exp) ** (-1./f)) ** (-f)  &
				/ solubility

	end function bubble_k
!
!
!
	real function sol_co2(sal, temp)
!===================================
! solubility of carbon dioxide in mol/(m^3 Pa)
! input: salinity in psu, temperature in K
! Weiss, R.F. (1974), Marine Chemistry, 2, 203-215
!
	real, parameter:: a1 = -58.0931, a2 = 90.5069, a3 = 22.2940
	real, parameter:: b1 = 0.027766, b2 = -0.025888, b3 = 0.0050578

	xlogk0 = a1 + a2*(100/temp) + a3 * log(temp/100) + sal * &
			(b1 + b2 * (temp/100) + b3 * (temp/100) *(temp/100))

! calc. of solubility and conversion from mol/(L atm) to mol/(m^3 Pa)
	sol_co2 = exp(xlogk0) / 101.325

	end function sol_co2

!
!
!
	real function diff_co2(temp)
!===============================
!...molecular gas diffusion coefficient for carbon dioxide in m^2/s
!...as function of temperature in K

	real, parameter:: c0 = 1.0e-9, c1 = 0.9e-9, c2 = 24., t0 = 273.15
	real temp
	diff_co2 = c0 + (temp - t0) * c1 / c2

	end function diff_co2
!
!
!
	real function foam(wind)
!===========================
! foam coverage according to Monahan and O'Muircheartaigh (1980)
! input: wind speed in m/s

	real wind
	real, parameter:: a1 = 3.84e-6, a2 = 3.41

	foam = a1 * wind ** a2

	end function foam
!
!
!
	real function vis_air(tair)
!==============================
! kinematic viscosity of air - Andreas (1989) crrel rep. 89-11, m^2/s
! input: temperature in K

	real tair
	real, parameter:: c1 = 1.326e-5, c2 = 6.542e-3, c3 = 8.301e-6, c4 = -4.84e-9
	real, parameter:: t0 = 273.15

	ta = tair - t0

	 vis_air = c1 * ( 1 + c2 * ta + c3 * ta * ta + c4 * ta * ta * ta)

	end function vis_air