$debug c Carol Anne's flux code sent from Judy Curry via Alam Afshan PROGRAM SHELL CALL FLUX(3,ARGV) END SUBROUTINE FLUX (ARGC,ARGV) C Arguments and Functions INTEGER ARGC CHARACTER ARGV REAL HUMLOW,getZoQ C ** Input variables REAL windsp,airT,qa,seaT,zu,zT,zQ,KV C** Other variables from TIWE data ** REAL jul,time,dir,qs,hsmall,hlc,tau,rs,rl,rain,rlup,htot, 1 direction C ** Output and other variables ** REAL uStar,thetas,qStar,zo,zL,Rr,Rt REAL Rq,skinT,timeout,nu,negHs,negHl,nuWat real virpottemp,ro,cp,lv,st,da,qnet,renewaltime C ** flux variables ** REAL M,Hs,Hl C ** Drag variables ** REAL Ri,uS,dU,dT,dQ,zoN,zoT,zoQ REAL PUZ,PTZ,PQZ,ZTL,ZQL REAL zLN,test1,test3 INTEGER n1,n3 C COMMON /FLUXM/Sc,Prt,KV,G,watTyp,BETAG,BETAC,BETAS,U_SFC,cmin, 1 hp,Prair,PrWat,nuWat,nu,CON_P,zotDet,zoDet,ID,BWU,zoinit,wa,w,zz, 1 wave_len C initialize DATA nu/1.5e-05/ DATA G/9.81000/ DATA KV/0.4/ DATA H/1000.0/ DATA B/0.019/ DATA SBOLTZ/5.67e-08/ DATA CON_P/0.70/ DATA Prt/0.85/ DATA Prair/0.7/ DATA PrWat/7.0/ DATA nuWat/1.5e-06/ DATA Sc/0.6/ DATA zoDet/5/ DATA zotDet/1/ DATA PI/3.14159/ data thermaldiff/2.e-5/ data moleculardiff/2.4e-5/ DATA watTyp/1/ DATA ID/0/ C INITIALIZE ZU,ZT,ZQ DATA zu/10./ DATA zT/11./ DATA zQ/11./ C c OPEN(UNIT=1,FILE='flux.astx.gmt',STATUS='old') c OPEN(UNIT=1,FILE='test.in',STATUS='old') c OPEN(UNIT=2,FILE='fluxout.two',STATUS='unknown') c c open input/output files and set up fixed instrument levels c open and read statements copied from coar2_5b.for c open(unit=3,file='d:\123\wave\test2_5b.dat') open(unit=4,file='d:\123\wave\test2_5b.out') c if(argc.NE.3) THEN STOP 10 END IF DO 100 J=1,999 c READ(1,*) time, uS, windSp, direction, seaT, airT, qs, qa, hs, c 1 hlc, tau, rs, rl, htot, rain c READ(3,505)ADUM ! skip header line of test data 505 FORMAT(A1) do while (.not. eof(3)) !start of loop index=index+1 !count data records (hours) c c read Chris Fairall's Moana Wave data from test file, 116 lines(hours) c hsb,hlb,tub are bulk fluxes calculated independently by Chris for comparison c use Hemantha's MSP data at 6m for Ts to demonstrate warm layer calculation c 700 READ(3,705,end=910)xtime,u,tt,t,q,hsb,hlb,tub, 2 rs,rl,rain,xlat,xlon,ts 705 FORMAT(f12.0,10F6.0,3F9.0) c c tt is Chris' floating temperature sensor at 0.05m depth ** Set Initial Variables ** uStar = 0.0 thetas = 0.0 qStar = 0.0 zo = 0.0005 zL = 1.0 Rr = 0.0 Rt = 0.0 Rq = 0.0 uS = 0.0 C ** Calculate the specific humidity of the water ** qs = humlow(seaT) C ** TIWE data has qa in g/kg: This program wants kg/kg ** qa = qa/1000 C ** Prepare some variables for iteration loop ** virPotTemp = ((airT + 273.16)*(1 + 0.608*qa)) - 273.16 ro = (3.48383e-03*1.01325e05)/(virPotTemp+273.16) Cp = 1.004*(1 + 0.9*qa)*1e03 Lv = 4.1868*(597.31-0.56525*airT)*1e03 c /** Prepare some variables for iteration loop **/ dU = windsp-uS dT = airT-seaT dQ = qa-qs uStar = 0.04*dU qStar = 0.003*dQ thetas = 0.003*dT rlup = 0.986*SBOLTZ*((seaT+273.15)**4.0) C** Begin iteration loop using Z/L as test ** DO 20 N3=1,50,1 C ** Begin iteration loop for zoN ** zo = getZo(uStar,windSp,qa,thetaS,airT,qStar,seaT,zu) PUZ = PSI(1,zL) ZTL = zL*zt/zu ZQL = zL*zq/zu PTZ = PSI(2,ZTL) PQZ = PSI(2,ZQL) uStar = dU*0.4/(log(zu/zo)-PUZ) Rr = zo*uStar/nu call getFlux(airT, qa, uStar, thetaS, qStar, M, Hs, Hl) Qnet = -rl+rlup + Hs + Hl renewalTime =getRenewalTime(zo,uStar,ro,Cp,1.0/(airT+273.15), 1 -Qnet, 1.4e-05) St = getStanton(uStar, renewalTime) Da = getDalton(uStar, renewalTime) zoT = getzoT(Rr, uStar, zo, St) zoQ = getzoQ(Rr, uStar, zo, Da) thetas = dT*0.4/(log(zt/zoT)-PTZ) qStar = dQ*0.4/(log(zq/zoQ)-PQZ) CALL getFlux(airT,qa,uStar,thetas,qStar,M,Hs,Hl) skinT =getSkin(zo,uStar,seaT,rl-rlup,-Hs,-Hl,rs) qs = humlow(skinT) dT = airT-skinT dQ = qa-qs rlup = 0.986*SBOLTZ*((skinT+273.15)**4.0) zLN = zeta(airT,qa,uStar,thetas,qStar,zu) test3 = ABS((zL - zLN)/(zL + 1.0e-08)) if(test3.LT.1.e-2) GO TO 220 zL = zLN IF(n3.EQ.50) THEN PRINT *,"Program fails to converge" PRINT *,"Program is exiting now" STOP 20 END IF 20 CONTINUE C** Now we compute final fluxes ** 220 CALL getFlux(airT,qa,uStar,thetas,qStar,M,Hs,Hl) print*,M, Hs, Hl, zL,zot,skinT-seaT,Rr READ (1,*,END=999) time, uS, windSp, direction, seaT, airT, qs, 1 qa, hs, hlc, tau, rs, rl, htot, rain 100 CONTINUE 999 PRINT *,'END OF DATA, PROGRAM TERMINATING' c CLOSE (2) CLOSE (1) STOP 999 END FUNCTION humlow(temp) REAL temp REAL wetTemp,A(8),x,es,qs,P INTEGER i,j C INITIALIZE DATA A(1)/4.436519e-01/ DATA A(2)/1.428946e-02/ DATA A(3)/2.650649e-04/ DATA A(4)/3.031240e-06/ DATA A(5)/2.034081e-08/ DATA A(6)/6.136821e-11/ DATA P/1008./ x=0. wetTemp = temp DO 30 i=1,6,1 j = 7-i x = (x+A(j))*wetTemp 30 CONTINUE es = (6.1078+x)*(1-2*0.0345*18/58.45) qs = 0.622*es/(P-0.378*es) humlow=qs RETURN end FUNCTION drag(u10) REAL u10 REAL CD,RAN(7),A(7),B1(7),P(7) INTEGER k,i COMMON /FLUXM/Sc,Prt,KV,G,watTyp,BETAG,BETAC,BETAS,U_SFC,cmin, 1 hp,Prair,PrWat,nuWat,nu,CON_P,zotDet,zoDet,ID,BWU,zoinit,wa,w,zz, 1 wave_len C C initialize values C DATA RAN(1)/2.2/ DATA RAN(2)/5.0/ DATA RAN(3)/8.0/ DATA RAN(4)/25.0/ DATA RAN(5)/50.0/ C DATA RAN(6)/50.0/ DATA A(1)/0.0/ DATA A(2)/0.771/ DATA A(3)/0.867/ DATA A(4)/1.2/ DATA A(5)/0.0/ DATA B1(1)/1.08/ DATA B1(2)/0.0858/ DATA B1(3)/0.0667/ DATA B1(4)/0.025/ DATA B1(5)/0.073/ DATA P(1)/-0.15/ DO 40 I=2,5 P(I) = 1.0 40 CONTINUE k = ID-2 if(k.LT.0 ) THEN if(u10.GT.50.0) THEN CD = 3.7e-03 elseif(u10.LT.0.3) THEN CD = 1.5e-03 else DO 120 I=1,7,1 IF(u10.LE.RAN(I)) then if (I.eq.1) then CD = (A(1)+B1(1)*((u10**P(1))))/1000. else GO TO 199 end if END IF CD = (A(I+1)+B1(I+1)*((u10**P(I+1))))/1000. 120 CONTINUE 199 CONTINUE END IF ELSEIF(k.EQ.0) THEN CD = (0.61+0.063*u10)/1000. else if(u10.GT.11.0) THEN CD = (0.49+0.065*u10)/1000. else CD = 1.2e-03 end if END IF drag=CD RETURN end FUNCTION PSI(iflag, zL) INTEGER iflag REAL psi,zL,chi REAL stabfn,F,psik,chic,chik,psic,A,B,C,D C initialize DATA A/1.0/ DATA B/0.667/ DATA C/6.0/ DATA D/0.35/ C/** This is basically original C version with some Fairall modification **/ C/* if(zL < 0.0) { C if(iflag != 1) { C chi = ((1.0 - 16.0*zL)**0.5) C stabfn = 2.0 * log((1.0 + chi)/2.0) C } C else { C chi = ((1.0 - 16.0*zL)**0.25) C stabfn = 2.0 * log((1.0 + chi)/2.0); C stabfn += log((1.0 + chi*chi)/2.0); C stabfn -= 2.0*atan(chi); C stabfn += PI/2.0; C } C** This is new Fairall version of convective form **/ if(zL.LT.0.0) THEN F = 1/(1 + zL*zL) chik = (1-16*zL)**0.25 if(iflag.NE.1) then psik = 2*log((1+chik*chik)/2.0) chic = (1-12.87*zL)**(1.0/3.0) psic = 1.5*log((chic*chic+chic+1.0)/3.0) psic = psic+(-(3.0**0.5)*atan((2*chic+1.0)/(3.**0.5))) psic = psic+(4*atan(1.0)/(3.**0.5)) stabfn = F*psik+(1-F)*psic else psik = 2*log(1+chik/2.0)+log((1+chik*chik)/2.0) psik = psik+(-2.0*atan(chik)+2*atan(1.0)) chic = (1-12.87*zL)**(1.0/3.0) psic = 1.5*log((chic*chic+chic+1.0)/3.0) psic =psic+(-(3.0**0.5)*atan((2*chic+1.0)/(3.0**0.5))) psic = psic+(4*atan(1.0)/(3.**0.5)) stabfn = F*psik+(1-F)*psic END IF else if(zL.EQ.0.0) THEN stabfn = 0.0 C* else { C stabfn = -6.0 * log(1.0 + zL) else if(iflag.EQ.1) THEN stabfn = -1.0*(A*zL+B*(zL-C/D)*exp(-D*zL)+B*C/D) else stabfn = (1+2*a*zL/3)**1.5+b*(zL-c/d)*exp(-d*zL) stabfn = stabfn+(b*c/d - 1) stabfn = -1.0*stabfn end if end if psi=stabfn RETURN end FUNCTION zeta(T,Q,USR,TSR,QSR,Z) REAL T,Q,USR,TSR,QSR,Z,von,g,ta,tv,tvsr,ob,zl,zeta c initialize DATA von/0.4/ DATA g/9.81/ ta = 273.16+T tv = ta*(1+0.61*Q) tvsr = TSR*(1.0+0.61*Q)+0.61*ta*QSR if(tvsr.EQ.0.0) THEN zl = 0.0 else ob = tv*USR*USR/(g*von*tvsr) zl = Z/ob END IF zeta=zl RETURN end FUNCTION getzoQ(Re,uStar,zo,dalton) REAL Re,uStar,zo,zoQ,dalton COMMON /FLUXM/Sc,Prt,KV,G,watTyp,BETAG,BETAC,BETAS,U_SFC,cmin, 1 hp,Prair,PrWat,nuWat,nu,CON_P,zotDet,zoDet,ID,BWU,zoinit,wa,w,zz, 1 wave_len zoQ = zo*exp(0.4*(5 - 1/(Prt*Dalton))) getzoQ=zoQ RETURN end FUNCTION getzoT(Re,uStar,zo,stanton) REAL Re,uStar,zo,zoT,stanton COMMON /FLUXM/Sc,Prt,KV,G,watTyp,BETAG,BETAC,BETAS,U_SFC,cmin, 1 hp,Prair,PrWat,nuWat,nu,CON_P,zotDet,zoDet,ID,BWU,zoinit,wa,w,zz, 1 wave_len zoT = zo*exp(0.4*(5 - 1/(Prt*Stanton))) getzoT=zoT RETURN end FUNCTION getZo(uStar,windSp,qa,thetaS,airT,qStar,seaT,zu) REAL uStar,windsp,qa,thetas,airT,qStar,seaT,zu,zo REAL u10,CD,C,zoN,test1,betag,betas,betac,U_SFC REAL wa, monin_inv, bStar, uStarOld, wStar2, dU2, test2 INTEGER n1 REAL KV C C function type declarations C REAL drag,getZoM,getbstar C COMMON /FLUXM/Sc,Prt,KV,G,watTyp,BETAG,BETAC,BETAS,U_SFC,cmin, 1 hp,Prair,PrWat,nuWat,nu,CON_P,zotDet,zoDet,ID,BWU,zoinit,wa,w,zz, 1 wave_len C C PRINT *,'ZODET =',zoDet if(zoDet.EQ.1.) THEN zo = 0.011*uStar*uStar/9.81 else if(zoDet.EQ.2.) THEN if(windsp.LT.10) THEN CD = 1.14E-3 else CD = 1.0E-3*(0.49 * 0.065*windsp) ENDIF zo = 10.0*exp(-0.4*(CD**(-0.5))) else if(zoDet.EQ.3.) THEN if(windsp.LT.2.2) THEN CD = 1.08E-3*(windsp**(-0.15)) else if(windsp.GE.2.2.AND.windsp.LT.5.0)THEN CD = (0.771 + 0.0858*windsp)*1.0E-3 else if(windsp.GE.5.0.AND.windsp.LT.8.0) THEN CD = (0.867 + 0.0667*windsp)*1.0E-3 else if(windsp.GE.8.0.AND.windsp.LT.25.0)THEN CD = (1.2 + 0.025*windsp)*1.0E-3 else CD = 0.073*windsp*1.0E-3 END IF zo = 10.0 * exp(-0.4*(CD**(-0.5))) else if(zoDet.EQ.4.) THEN C PRINT *,'REACHED START OF LOOP 60' zo=0.0005 c PRINT *,' zo=',zo DO 60 n1=1,50 u10 = uStar*log(10.0/zo)/0.4 CD = drag(u10) C = 1.0/sqrt(CD) zoN = 10.0/exp(0.4*C) test1 = ABS((zoN - zo)/(zo+ 1.0e-08)) C print *,'test1=',test1 if(test1.LT.0.01) GO TO 61 zo= zoN if(n1.EQ.50) then PRINT *,"No convergence of CD after 50 runs" PRINT *,"Program exiting..." STOP60 END IF 60 CONTINUE 61 continue C zo=zoN C print *, 'reached end of loop 60' else if(zoDet.EQ.5.) THEN betag = 1.0 betac = 1.0 betas = 0.0 wa = 20.0 u_sfc = 0.55*uStar monin_inv = 0.0 n1 = 0 bStar = getbStar(qa, thetaS, airT, qStar) monin_inv = KV*bStar/(uStar*uStar) if(abs(monin_inv).lt.0.001) monin_inv = 0.001*monin_inv/ 1 abs(monin_inv) if(abs(monin_inv).gt.10.00) monin_inv = 10.00*monin_inv/ 1 abs(monin_inv) zo = getZoM(seaT,airt,windSp,zu,uStar,monin_inv) uStarOld = uStar wStar2 = (abs(uStar*bStar)*1000.0)**0.666667 dU2 = (windSp - u_sfc)*(windSp - u_sfc)+CON_P*wStar2 if(dU2.lt.0.000001) dU2 = 0.000001 uStar = sqrt(dU2)*0.4/(log(zu/zo + 1)) test2 = abs((uStar - uStarOld)/uStar) c print*,zo,seat,windsp,zu,ustar,'zo:after getzom' end if getZo=zo RETURN end SUBROUTINE getFlux(airT,qa,uStar,thetas,qStar,M,Hs,Hl) C This subroutine computes the surface momentum, heat and moisture fluxes REAL airT,qa,uStar,thetas,qStar,M,Hs,Hl REAL virPotTemp,ro,Cp,Lv virPotTemp = ((airT+273.16)*(1+0.608*qa))-273.16 ro = (3.48383e-03*1.01325e05)/(virPotTemp+273.16) Cp = 1.004*(1+0.9*qa)*1e03 Lv = 4.1868*(597.31-0.56525*airT)*1e03 M = -1*ro*uStar*uStar Hs = -1*ro*Cp*thetas*uStar Hl = -1*ro*qStar*uStar*Lv RETURN end FUNCTION getSolar(solar,watTyp) REAL solar,gamma1,gamma2,R,z,solarAbs if(watTyp.EQ.0)THEN getsolar = 0.0 RETURN else if(watTyp.EQ.1) THEN R = 0.58 gamma1 = 0.35 gamma2 = 23 else if(watTyp.EQ.1.1) THEN R = 0.62 gamma1 = 0.60 gamma2 = 20 else if(watTyp.EQ.1.2) THEN R = 0.67 gamma1 = 1.00 gamma2 = 17 else if(watTyp.EQ.2) then R = 0.77 gamma1 = 1.50 gamma2 = 14 else if(watTyp.EQ.3) THEN R = 0.78 gamma1 = 1.40 gamma2 = 7.9 END IF z = 0.001 solarAbs = solar-solar*(R*exp(-z/gamma1)+(1-R)*exp(-z/gamma2)) getSolar=solarAbs RETURN END FUNCTION getbStar(qa, thetaS, airT, qStar) REAL qa, thetas, airt, qstar, bstar BSTAR = 9.81 * (( 1.0000000 + 0.6100000 * QA ) * Thetas + & 0.61000000 * (Airt + 273.15) * QSTAR ) / & ( Airt + 273.15 + 0.61000 * QA ) if ( abs( bstar ).lt.0.000001 ) then if ( bstar.eq.0.0 ) then bstar = 0.000001 else bstar = 0.000001 * bstar / abs( bstar ) endif else endif getbstar =bstar RETURN end FUNCTION PHIU(zref,MONINV,zz) REAL zref,MONINV,zz,PHIU,zeta,zeta0 if(MONINV.LT.0) THEN zeta = (1.0 - 15.0*zref*MONINV)**0.25 zeta0 = (1.0 - 15.0*zref/zz*MONINV)**0.25 a1=zeta0*zeta0+1.0 a2=(zeta0+1.0)**2. b1=zeta*zeta+1.0 b2=(zeta+1.0)**2. PHIU = log((a1*a2)/(b1*b2)) 2 +2.0*(atan(zeta)-atan(zeta0)) else PHIU = 7.0*zref*MONINV END IF RETURN end FUNCTION getZoM(Tw,ta,u,zref,ustar,moninv) LOGICAL DONE real B, BETAC, BETAG, BETAS, CRIT, & DENWAT, G, H, KV, MS, MW, NU, PI2, S, & SFCTEN, TA, TW, U, U_SFC, USTAR, ZREF, WA, ZZ, W,w_eql, & WA_OLD, hp,hp_old,period,cmin,moninv,dum,hzzg,phiu,z0 COMMON /FLUXM/Sc,Prt,KV,G,watTyp,BETAG,BETAC,BETAS,U_SFC,cmin, 1 hp,Prair,PrWat,nuWat,nu,CON_P,zotDet,zoDet,ID,BWU,zoinit,wa,w,zz, 1 wave_len w = 1.000 w_eql = 0.81 zz = 10000.0 s=.0349 C convergence critereon (fractional change) CRIT = 1.e-3 C equilibrium constant, b, for capillary waves [] B = 0.01900000 C von Karman's constant for momentum [] KV = 0.40 C gravitational acceleration (at sea level) [m/s^2] G = 9.8100000 C approximate height of boundary layer (at low wind speed) H = 1000.0 C molecular mass of salt [kg/kmol] MS = 58.443000 C molecular mass of water [kg/kmol] MW = 18.016000 C two PI PI2 = 6.28318 C determine the viscosity of air c NU = 1.3260000e-5 * (1.00000 + 0.006542 * TA + 8.301000E-6 c & * TA * TA + 4.840000E-9 * TA * TA * TA ) C determine the surface tension SFCTEN = 7.6100000E-2 - 1.55000E-4 * TW + 2.77000E-5 * S * MS / MW C determine the density of pure water DENWAT = (999.8396 + 18.224944 * TW - 0.007922210 * TW * TW ) / & ( 1.0000000 + 0.018159725 * TW ) C determine the density of saline water VA = ( 12.97000 + 0.23400000 * TW - 4.2100000E-3 * TW * TW + & 2.8570000E-5 * TW *TW *TW ) * 10.00000E-3 + & SQRT( S * DENWAT * 10.00000E-3 ) * & 2.9820000E-3 - 4.9700000E-5 * TW + 6.0320000E-7 * TW * TW DENWAT = DENWAT * ( 1.000 + S ) / ( 1.000 + VA * DENWAT * S / MS ) U_SFC = 0.55 * USTAR denwat = 1027.0 sfcten = (-1.55e-04*(273.15+Tw)) + 0.118 CMIN = SQRT( 2.0 * SQRT( G * SFCTEN / DENWAT ) ) 410 continue C determine if the phase speed of the dominant waves is physically C reasonable. If the phase speed is less than the physical minimum, C the assume that assume the surface is smooth. IF ( BETAG .NE. 0.0 .AND. ( WA * USTAR .LT. CMIN ) ) THEN BETAG = 0.0 BETAC = 0.0 BETAS = 1.0 ENDIF CMIN = SQRT( 2.0 * SQRT( G * SFCTEN / DENWAT ) ) C 'Guess' at roughness length Z0 = BETAS * 0.1100000 * NU / abs(USTAR) + & BETAC * B * SFCTEN / ( USTAR * USTAR * DENWAT ) + & 0.480 / WA * USTAR * USTAR / G C determine the wave characteristics dum = ustar * ustar * wa * wa if(dum.gt.cmin*cmin) then wave_len = 3.14159 * (dum+(sqrt(dum * dum - cmin**4.0))) / G else wave_len = 3.14159*(2*dum)/G endif period = wave_len / (wa * ustar) hp =0.062* SQRT( G * USTAR * PERIOD ** 3 ) c /* limit the height of breaking waves */ IF ( HP .GT. 0.137 * WAVE_LEN ) HP = 0.137 * WAVE_LEN c /* determine the mean water surface speed */ u_sfc = betag * ( 0.5 * PI2 * hp / wave_len )**2. * ustar * wa + 1 ( 1.0 - betag ) * 0.55 * ustar c /* determine wave age */ clayson should use the following for phi_u as function phiu requires 3 arguments c phi_u = phiu(betag * hp , moninv, hp / z0 ) phi_u = phiu(betag * hp , hp / z0 ) hzzg = hp * zz / zref if ( hzzg.lt.0.00001 ) hzzg = 0.00001 wa = (log(hzzg + 1.0) + phi_u ) / (0.4 * w * w_eql) c /* insure that Wa is not horribly unreasonable */ if ( wa.gt.150.0 ) wa = 150.0 if ( wa.lt.1.1 ) wa = 1.08 c /* determine zref / z0 */ zz = zref / ( ( betas * 0.1100 * nu / abs( ustar )) + 1 betac * B * sfcten / ( ustar * ustar * denwat ) + betag * 1 0.480 / wa * ustar * ustar / G ) c /* insure that the roughness length is not too absurd */ if ( zref.gt.zz * 25 ) zz = zref / 25.0 if( zref/zz.gt.1e-03) zz = zref/1e-03 if ( betag.ne.0.and.(wa * ustar.lt.cmin)) go to 410 getZoM =zref/zz RETURN C end of subroutine z_O_zo END function getSkin(zo, uStar, bulkTemp, rl, Hs, Hl,rs) real zo,uStar,bulkTemp,rl,Hs,Hl,rs real Qsolar,Qnet,roWater,CpWater,uStarWater,salinity, 1 roAir,skinTemperature,kappawater,thermalExpansion,viscosity, 1 renewalTime salinity = 35 roWater = 1021.7 CpWater = 4002.0 roAir = 1.25 kappawater = 1.4e-07 thermalExpansion = 3196e-07 viscosity = 1e-06 c /** get solar radiation absorbed over skin **/ Qsolar = getSolar(rs,wattyp) c /** Get Net Heat Loss from ocean **/ Qnet = rl + Hs + Hl + Qsolar c /** get u* for water **/ uStarWater = sqrt(roAir/roWater)*uStar c /** now get skin temperature **/ renewalTime = getRenewalTime(zo,uStarWater, roWater, CpWater, 1 thermalExpansion, Qnet, viscosity) skinTemperature = sqrt(renewalTime)*Qnet/(roWater*CpWater* 1 sqrt(kappawater)) getskin = skinTemperature + bulkTemp return end function getRenewalTime(zo, uStar, density, cp, alpha, Qnet, 1 viscosity) real zo, uStar, density, cp, alpha, Qnet, viscosity real renewalTime,Rfcr,Rfo,Cshear,Cconv,shearContribution, 1 convContribution Rfcr = -2.0e-04 if(density.lt.50) then Rfcr = -2.0e-01 c changed to 53.5824 = 7.32**2 c Cshear = 53 Cshear = 53.5824 Cconv = 0.80 else Cshear = 209 Cconv = 3.13 endif Rfo=alpha*9.81*Qnet*viscosity/(density*cp*uStar**4.0) if(Rfo.eq.abs(Rfo)) Rfo = rfo*-1.0 c /** now get time rate **/ shearContribution = Cshear*sqrt(viscosity*zo/(uStar**3.0)) convContribution = 1 Cconv*sqrt(viscosity*density*cp/(alpha*9.81*abs(Qnet))) renewalTime = (convContribution-shearcontribution)*exp(-Rfcr/Rfo) renewalTime = renewaltime + shearContribution getrenewalTime = renewaltime return end function getStanton(uStar, renewalTime) real uStar,renewalTime,getStanton thermaldiff = 2.e-5 getStanton = sqrt(thermalDiff/(uStar*uStar*renewalTime)) return end function getDalton(uStar, renewalTime) real uStar,renewalTime,getDalton,moleculardiff moleculardiff = 2.4e-5 getDalton = sqrt(molecularDiff/(uStar*uStar*renewalTime)) return end