SUBROUTINE BULK(SWRAD,EMP,RELAX) C***********************************************************************C C SUBROUTINE TO CALCULATE SURFACE FLUXES FROM ETA ATMOS. MODEL DATA C C WUSURF, WVSURF, WTSURF (LOSS), EMP ARE CALCULATED IN BULK C C SWRAD IS GIVEN DIRECTLY FROM ETA MODEL C C DOWNWARD LONGWAVE RADIATION IS GIVEN DIRECTLY FROM ETA MODEL C C UPWARD LONGWAVE RADIATION IS GIVEN BY BIGNAMI 1995 C UNITS IN MKS C C FORMULAS TAKEN FROM CASTELARI ET.AL. 1992 C C***********************************************************************C c-------------------------------------------------------------------- c coefficients ( in MKS ) : c----------------------------------------------------------------- c c --- surface air pressure, expsi, dry air gas constant c data ps,expsi,rd / 1013., 0.622, 287./ c c --- turbulent exchange coefficients ( from Rosati & Miyakoda 1988 ) c data ce1,ch1 / 1.1e-3, 1.1e-3/ c c --- turbulent exchange coefficients ( from Budyko 1963 ) c data ce2,ch2 / 2.1e-3, 2.1e-3/ c c --- air density, Stefan-Boltzmann constant , ocean emissivity c data arho,sigma ,emiss /1.2, 5.67e-8, .97/ c c --- Solar constant , specific heat capacity c data solar ,cp /1350., 1005./ c c data ckelv /273.16/ c const = 0.622/1013. pi = acos(-1.) c c---------------------------------------------- sqe = 0.e0 sqh = 0.e0 ssol = 0.e0 sqbw = 0.e0 sqbwup = 0.e0 swind = 0.e0 totcell = 0.e0 stemp = 0.e0 srain= 0.e0 sevap= 0.e0 ssstk= 0.e0 stnowk= 0.e0 srhum = 0.e0 c cd_max = -1.e6 cd_min = 1.e6 t_max = -1.e6 t_min = 1.e6 r_max = -1.e6 r_min = 1.e6 tair_max = -1.e6 tair_min = 1.e6 tsst_max = -1.e6 tsst_min = 1.e6 c-------- -------- -------- -------- -------- -------- -------- -------- c DO 1000 i=1,im DO 1000 j=1,jm if (fsm(i,j).ne.0.) then unow = uair(i,j) vnow = vair(i,j) tnow = tair(i,j) rhnow = rhum(i,j) precip= rain(i,j) sst_o = t(i,j,1)+10. dlong = rlond(i,j) sola = sol(i,j) c c c c --- compute sp for bulk coeff. c SP = sqrt(unow*unow+vnow*vnow) c c --- SST data converted in Kelvin degrees c --- TNOW is already in Kelvin c sstk = sst_o + ckelv tnowk = tnow c c c ---calculates the Saturation Vapor Pressure at air temp. and at sea temp. c ---esat(Ta) , esat(Ts) c esatair = esk(tnowk) esatoce = esk(sstk) c c --- calculates the saturation mixing ratios at air temp. and sea temp. c --- wsat(Ta) , wsat(Ts) c wsatair = (expsi/ps) * esatair wsatoce = (expsi/ps) * esatoce c c --- calculates the mixing ratio of the air c --- w(Ta) c wair = 0.01 * rhnow * wsatair c c --- calculates the density of the moist air c rhom = 100.*(ps/rd) * (expsi*(1.+wair)/(tnowk*(expsi+wair))) c c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c QBW (the long wave flux) c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c c ----- Bignami et al. 1995 QBW=0.98*sigma*sstk**4 - dlong c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c --- calculates the term : ( Ts - Ta ) deltemp = sstk - tnowk c c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c QH (the sensible heat flux) and QE (the latent heat flux) c---- ----- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c c ---- Kondo scheme ss=0. if(sp.gt.0.0)ss=deltemp/(sp**2.) c c --- calculate the Stability Parameter : c stp=ss*abs(ss)/(abs(ss)+0.01) c c --- for stable condition : IF(ss.lt.0.) THEN if((stp.gt.-3.3).and.(stp.lt.0.)) then fh=0.1+0.03*stp+0.9*exp(4.8*stp) fe=fh else if(stp.le.-3.3) then fh=0. fe=fh c endif endif c c --- for unstable condition : ELSE fh=1.0+0.63*sqrt(stp) fe=fh ENDIF c --- calculate the coefficient CH3,CE3,CD3 c ch2=1.3e-03*fh ce2=1.5e-03*fe c--------------------------------------------- c--------------------------------------------- c --- calculates the QH and QE with the constant values for Ch,Ce from c --- Budyko (1963) c QH = rhom*cp*ch2*sp*deltemp c c --- calculates the term : esat(Ts)-r*esat(Ta) c evap = esatoce - rhnow*0.01*esatair c c --- calculates the term : Ce*|V|*[esat(Ts)-r*esat(Ta)]0.622/1013 c --- Evaporation rate [kg/(m2*sec)] c EVAP = rhom*ce2*sp*evap*const c c --- calculates the LATENT HEAT FLUX QE in MKS ( watt/m*m ) c --- QE = L*rho*Ce*|V|*[esat(Ts)-r*esat(Ta)]0.622/1013 c QE = evap*heatlat(sst_o) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c WFLUX ( the water flux ) in m/sec c---- -- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c WFLUX = evap/1023. - precip c c ---- here relax surface salinity to monthly climatology c --- Reverse Sign for model needs relax(i,j) = 0.0 emp(i,j)= -wflux c c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c QU ( the net upward flux ) c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c c --- calculates : Qu = Qb + QH + QE ( Rosati,Miyakoda 1988 ; eq. 3.2 ) c QU = qbw + qh + qe c c --- Devide by 1023*3991 (rho*cp) for model needs c --- DO NOT Reverse Sign because it is already positive c wtsurf(i,j) = QU*(1./4.082793e6) swrad(i,j) = -sola/4.082793e6 c c---- ----- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c TAUX, TAUY ( the wind stresses ) c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c c --- calculates the Drag Coefficient Cd with variable value from c --- Hellerman & Rosenstein (1983) c c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ c --- (5-10-92) we discover that CD of H. R. needs to be used with c --- deltem = Ts - Ta c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ cd1 = cd(sp,deltemp) if (cd1.gt.1.3e-3) cd1=1.3e-3 c c --- calculates the wind stresses in MKS ( newton/m*m ) c --- taux= rho*Cd*|V|u tauy= rho*Cd*|V|v c TAUX= rhom*cd1*sp*unow TAUY= rhom*cd1*sp*vnow c --- Reverse Sign and devide by 1023 (rho) for model needs wusurf(i,j)=-taux/1023. wvsurf(i,j)=-tauy/1023. c c ------- MULTIPLY BY MASK ------------------------------------- c wusurf(i,j)=wusurf(i,j)*fsm(i,j) wvsurf(i,j)=wvsurf(i,j)*fsm(i,j) wtsurf(i,j)=wtsurf(i,j)*fsm(i,j) emp(i,j)=emp(i,j)*fsm(i,j) cell = fsm(i,j)*dx(i,j)*dy(i,j) totcell = totcell + cell sqbw = sqbw + qbw*cell ssol = ssol + sola*cell sqe = sqe + qe*cell sqh = sqh + qh*cell ssstk = ssstk + sstk*cell stnowk = stnowk + tnowk*cell swind = swind + cell*sqrt(unow*unow + vnow*vnow) srain = srain + cell*precip srhum = srhum + cell*rhnow*0.01 sevap = sevap + cell*evap/1023. c end if 1000 CONTINUE c ska=0. c do i=1,im c do j=1,jm c skat(i,j)=((-wusurf(i,j)*1023.)**2+ c .(-wvsurf(i,j)*1023.)**2)**0.5 c ska=max(ska,skat(i,j)) c enddo c enddo c do i=1,im c do j=1,jm c if (ska.eq.skat(i,j)) then c write(65,*) ska,-wusurf(i,j)*1023., c .-wvsurf(i,j)*1023.,i,j,time c endif c enddo c enddo c sqbw = sqbw/totcell ssol = ssol/totcell sqe = sqe/totcell sqh = sqh/totcell swind = swind/totcell srhum = srhum/totcell ssstk = ssstk/totcell stnowk = stnowk/totcell srain = srain/totcell *24.*3600.*360. sevap = sevap/totcell *24.*3600.*360. write(*,'(a4,11f7.2)') 'FLUX',time2,ssol,sqe,sqbw,sqh, 1 stnowk,ssstk,srhum,sevap,srain,swind c c -------------------------------------------------- RETURN END c function HEATLAT(t) c c --- calculates the Latent Heat of Vaporization ( J/kg ) as function of c --- the temperature ( Celsius degrees ) c --- ( from A. Gill pag. 607 ) c c --- Constant Latent Heat of Vaporization ( Rosati,Miyakoda 1988 ) c L = 2.501e+6 (MKS) c heatlat = 2.5008e+6 -2.3e+3*t c return end c function CD(sp,delt) c c --- calculates the Drag Coefficient as a function of the abs. value of the c --- wind velocity c --- ( Hellermann and Rosenstein ) c dimension a(6) data a / 0.934e-3,0.788e-4,0.868e-4 & ,-0.616e-6,-.120e-5,-.214e-5/ c cd = a(1) + a(2)*sp + a(3)*delt + a(4)*sp*sp & + a(5)*delt*delt + a(6)*sp*delt c return end c c============================================================================== c real function ESK(t) c c --- compute the saturation water vapor pressure from c --- temperature in kelvin (Lowe,1977; JAM,16,100,1977) c dimension a(7) data a /6984.505294,-188.9039310,2.133357675,-1.288580973e-2, & 4.393587233e-5,-8.023923082e-8,6.136820929e-11 / c esk = a(1) +t*(a(2)+t*(a(3)+t*(a(4)+t*(a(5)+t*(a(6)+t*a(7)))))) c return end c function CD2(sp,delt) dimension a(3),b(3) data a/0.862,0.088,0.00089/ data b/0.1034,0.00678,0.0001147/ w = max(2.5,min(32.5,sp)) cd0 = 1.e-3*(a(1)+a(2)*W-a(3)*W*W) cd1=1.e-3*(b(1)-b(2)*W+b(3)*W*W) CD2= CD0 + CD1*delt return end