C C SUBROUTINE BULK(IM,JM,TBIAS,FSM,TSURF,ALON,ALAT,IYR,IMO,IDAY, 1 IHOUR,IMIN,WUSURF,WVSURF,WTSURF,SWRAD,PME,UAIR,VAIR,TAIR, 2 RHUM,RAIN,CLOUD) C******************************************************************************* C THIS SUBROUTINE PROVIDES SURFACE BOUNDARY CONDITIONS FOR MOMENTUM, C HEAT AND SALT EQUATIONS SOLVED BY THE HYDRODYNAMIC MODEL IN CASES WHEN C THE ATMOSPHERIC MODEL PROVIDES CLOUD COVER DATA INSTEAD OF THE NET SOLAR C RADIATION FLUX AND THE DOWNWARD LONGWAVE RADIATION FLUX AS IT IS ASSUMED C IN THE FIRST VERSION OF BULK CODE. THE NET SOLAR RADIATION IS CALCULATED C ACCORDING TO THE REED FORMULA WHILE THE NET LONGWAVE RADIATION CAN BE C CALCULATED ACCORDING TO BIGNAMI OR MAY FORMULA (SEE LOGICAL VARIABLES C BIGNAMI_FORMULA & MAY_FORMULA BELOW) C C MOMENTUM, HEAT AND FRESHWATER FLUXES ARE CALCULATED FROM THE ATMOSPHERIC C PARAMETERS (WIND VELOCITY, AIR TEMPERATURE, RELATIVE HUMIDITY, PRECIPITATION, C CLOUD COVERAGE) PROVIDED BY THE C WEATHER PREDICTION MODEL AND THE MODEL'S SST USING PROPER AIR-SEA BULK C FORMULAE (Castellari et al., 1998, Korres and Lascaratos 2003). C ( Castellari et al., 1998. Journal of Marine Systems, 18, 89-114 ; C Korres and Lascaratos, 2003. Annales Geophysicae, 21, 205-220.) C C ALL UNITS ARE S.I. (M.K.S.) C C THE USER HAS THE OPTION TO CALCULATE THE NET LONGWAVE RADIATION FLUX C ACCORDING TO BIGNAMI OR MAY FORMULA. THIS IS DONE THROUGH THE LOGICAL C VARIABLES BIGNAMI_FORMULA & MAY_FORMULA C C SUBROUTINE BULK PROVIDES ITS OUTPUT INTO THE ARRAYS: C C 1. WUSURF(): X-COMPONENT OF WIND STRESS DIVIDED BY (-1)*RHO C 2. WVSURF(): Y-COMPONENT OF WIND STRESS DIVIDED BY (-1)*RHO C 3. SWRAD() : NET SOLAR RADIATION HEAT FLUX DIVIDED BY (-1)*RHO*Cpw C 4. WTSURF(): NET UPWARD HEAT FLUX DIVIDED BY RHO*Cpw C 5. PME() : PRECIPITATION MINUS EVAPORATION RATE C ( RHO: sea water density, Cpw: specific heat capacity of seawater ) C C C SUBROUTINE BULK NEEDS THE FOLLOWING INPUT: C C C MODEL RELATED DATA C 1. IM : NUMBER OF GRID POINTS IN X C 2. JM : NUMBER OF GRID POINTS IN Y C 3. TBIAS : CONSTANT VALUE SUBTRACTED FROM MODEL'S TEMPERATURE FIELD C 4. FSM() : THE MODEL MASK (1.:SEA GRID POINT, 0.:LAND GRID POINT) C 5. TSURF(): MODEL'S TEMPERATURE FIELD AT THE TOP VERTICAL LEVEL AND AT THE C CENTRAL TIME LEVEL C 6. ALON() : LOGNITUDE OF GRID POINTS C 7. ALAT() : LATITUDE OF GRID POINTS C 9. IYR : INTEGER VALUE CORRESPONDING TO CURRENT YEAR(for example iyr=2002) C 10. IMO : INTEGER VALUE CORRESPONDING TO CURRENT MONTH (1-->12) C 11. IDAY : INTEGER VALUE CORRESPONDING TO CURRENT DAY (1-->31) C 12. IHOUR : INTEGER VALUE CORRESPONDING TO CURRENT HOUR (0-->23) C 13: IMIN : INEGER VALUE CORRESPONDING TO CURRENT MINUTES (0--59) C C ATMOSPHERIC DATA: C 6. UAIR() : X-COMPONENT OF AIR VELOCITY (in m/s) AT 10m ABOVE SEA SURFACE C 7. VAIR() : Y-COMPONENT OF AIR VELOCITY (in m/s) AT 10m ABOVE SEA SURFACE C 8. TAIR() : AIR TEMPERATURE (in deg Kelvin) AT 2m ABOVE SEA SURFACE C 9. RHUM() : RELATIVE HUMIDITY (%) AT 2m ABOVE SEA SURFACE C 10. RAIN() : PRECIPITATION RATE (in m/s) C 11. CLOUD() : CLOUD COVERAGE IN TENTHS (0.-->1.) C C C Important: C SUBROUTINE BULK REQUIRES THAT THE ATMOSPHERIC DATA ARE ALREADY C INTERPOLATED IN TIME AND SPACE (i.e. mapped onto the model grid and C interpolated to the current time step. The user has to write his/her own C subroutines in order to map the raw atmospheric data onto the model grid C and interpolate them to the current time step) C C******************************************************************************* LOGICAL BIGNAMI_FORMULA, MAY_FORMULA PARAMETER (BIGNAMI_FORMULA=.TRUE., MAY_FORMULA=.FALSE.) DIMENSION FSM(IM,JM),PME(IM,JM),SWRAD(IM,JM),WUSURF(IM,JM), 1 WVSURF(IM,JM),WTSURF(IM,JM),TSURF(IM,JM),ALON(IM,JM), 2 ALAT(IM,JM) DIMENSION UAIR(IM,JM),VAIR(IM,JM),TAIR(IM,JM),RHUM(IM,JM), 1 RAIN(IM,JM),CLOUD(IM,JM) c c-------------------------------------------------------------------- c coefficients ( in MKS ) : c----------------------------------------------------------------- c c --- Sea water density data rho/1023./ c --- Sea water density times the specific heat of seawater data rho_cpw/4.082793e6/ c --- surface air pressure, expsi, dry air gas constant c data ps,expsi,rd / 1013., 0.622, 287./ 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 of air c data solar ,cp /1350., 1005./ c c data ckelv /273.16/ c const = 0.622/1013. 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) cld = cloud(i,j) sst_model = tsurf(i,j)+tbias C SST_MODEL IS THE MODEL'S TEMPERATURE (in deg Celsius) AT THE TOP LEVEL c c --- compute wind speed magnitude 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_model + 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 moist air c rhom = 100.*(ps/rd) * (expsi*(1.+wair)/(tnowk*(expsi+wair))) c c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c Calculate the net longwave radiation flux at the sea surface (QBW) c according to Bignami (Bignami et al., 1995) or May formula (May,1986) c c Bignami et al., 1995: Longwave radiation budget in the Mediterranean c Sea, J.Geophys Res., 100, 2501-2514. c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c ea12 = 0.01*rhnow*esatair if (bignami_formula) then QBW=0.98*sigma*sstk**4 -sigma*tnowk**4*(0.653+0.00535*ea12) @ *(1.+0.1762*cld*cld) end if if (may_formula) then QBW = (1.-0.75*(cld**3.4)) & * (sigma*(tnowk**4.)*(0.4 -0.05*sqrt(ea12)) & + 4.*sigma*(tnowk**3.)*(sstk-tnowk)) end if c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c --- calculate the term : ( Ts - Ta ) deltemp = sstk - tnowk c c Calculate turbulent exchange coefficients according to Kondo scheme c ( Kondo, J., 1975. Boundary Layer Meteorology, 9, 91-112) 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 ch2=1.3e-03*fh ce2=1.5e-03*fe c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c Calculate the sensible (QH) latent (QE) heat flux c and the evaporation rate (EVAP) c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 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 --- calculate 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 --- calculate the LATENT HEAT FLUX QE in MKS ( watt/m*m ) c --- QE = L*rhom*Ce*|V|*[esat(Ts)-r*esat(Ta)]0.622/1013 c QE = evap*heatlat(sst_model) c c c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c Calculate the water flux (WFLUX) in m/sec c---- -- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c WFLUX = evap/rho - precip pme(i,j)= -wflux c C Important note for Princeton Ocean Model users: C THE SALT FLUX ( WSSURF() ) IN POM MAIN CODE (REQUIRED FOR PROFT) SHOULD C BE CALCULATED AS: C do 3072 j = 1, jm C do 3072 i = 1, im C 3072 WSSURF(I,J)=PME(I,J)*(VF(I,J,1)+SBIAS) c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c Calculate the net upward flux (QU) at the sea surface c---- --- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c c --- calculates : Qu = Qb + QH + QE c QU = qbw + qh + qe c c --- 1. Devide upward heat flux by rho*Cpw c wtsurf(i,j) = QU/rho_cpw c --- Calculate net solar radiation flux according to Reed (Reed,1977) formula c c Reed, R.K.,1977: On estimating insolation over the ocean, J.Phys. c Oceanogr. 17, 854-871. alonp = alon(i,j) alatp = alat(i,j) call sol_rad(sol_net,cld,alonp,alatp,iyr,imo,iday,ihour,imin) c --- 1. Devide net solar radiation flux by rho*Cpw and reverse sign swrad(i,j) = -sol_net/rho_cpw c c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c Calculate the wind stress components (TAUX, TAUY) c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c c --- Calculate the Drag Coefficient Cd accoriding to c Hellerman & Rosenstein (1983) c cd1 = cd(sp,deltemp) c --- Calculate the wind stresses in MKS ( newton/m*m ) c --- taux= rhom*Cd*|V|u tauy= rhom*Cd*|V|v c TAUX= rhom*cd1*sp*unow TAUY= rhom*cd1*sp*vnow c --- Reverse Sign and devide by sea water density wusurf(i,j)=-taux/rho wvsurf(i,j)=-tauy/rho c end if c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- c Multiply all fluxes by model mask c---- ------ ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 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) pme(i,j)=pme(i,j)*fsm(i,j) 1000 CONTINUE c 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 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, 1983. JPO, 13,1093-1104) 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 SUBROUTINE SOL_RAD(SOLA,CLD,ALONP,ALATP,IYR,IMT,IDY,IHR,IME) C-------------------------------------------------------------------- DATA PI/3.141593/,RAD/.01745329/ c c c Panos scale solar radiation scal=1.0 c c --- Rosati and Miyakoda : c blon = alonp*pi/180. blat = alatp*pi/180. acl=cld qsw=qshort(iyr,imt,idy,ihr,ime,blat,blon,acl) sola = qsw c return end function qshort(iyr,imt,idy,ihr,ime,alat,alon,acl) c parameter(pi=3.1415927,degrad=pi/180.,raddeg=180./pi, * eclips=23.439*degrad) c dimension alpham(12),alb1(20),za(20),dza(19) c c --- alat,alon - (lat, lon) in radians !! c data solar/1350./ data tau /0.7/ data aozone /0.09/ data yrdays /365./ data alb1/.719, .656, .603, .480, .385, .300, .250, .193, .164 $ ,.131 , .103, .084, .071, .061, .054, .039, .036, .032, .031 $ ,.030 / c data za/ 90., 88., 86., 84., 82., 80., 78., 76., 74., 70., 66. $ ,62., 58., 54., 50., 40., 30., 20., 10., 0.0 / c data dza/8*2.0, 6*4.0, 5*10.0/ c c --- albedo monthly values from Payne (1972) as means of the values c --- at 40N and 30N for the Atlantic Ocean ( hence the same latitudinal c --- band of the Mediterranean Sea ) : c data alpham /0.09,0.08,0.06,0.06,0.06,0.06,0.06,0.06, & 0.06,0.07,0.09,0.10/ c c--------------------- calculations start ----------------------------- c c --- sun hour angle : c DTOR=DEGRAD XLCT=(1.0*IHR)+(1.0*IME/60.0) UT=XLCT IF (IMT.GT.2) THEN IYR1=IYR IMT1=IMT-3 ELSE IYR1=IYR-1 IMT1=IMT+9 ENDIF INTT1=INT(30.6*IMT1+0.5) INTT2=INT(365.25*(IYR1-1976)) SMLT=((UT/24.0)+IDY+INTT1+INTT2-8707.5)/36525.0 EPSILN=23.4393-0.013*SMLT CAPG=357.528+35999.050*SMLT IF(CAPG.GT.360.0) THEN G360=CAPG-INT(CAPG/360.0)*360.0 ELSE G360=CAPG ENDIF CAPC=1.915*SIN(G360*DTOR)+0.020*SIN(2*G360*DTOR) CAPL=280.460+36000.770*SMLT+CAPC IF(CAPL.GT.360.0) THEN XL360=CAPL-INT(CAPL/360.0)*360.0 ELSE XL360=CAPL ENDIF ALPHA=XL360-2.466*SIN(2*XL360*DTOR)+0.053*SIN(4*XL360*DTOR) GHA=15.0*UT-180.0-CAPC+XL360-ALPHA IF(GHA.GT.360.0) THEN GHA360=GHA-INT(GHA/360.0)*360.0 ELSE GHA360=GHA ENDIF DEC=ATAN(TAN(EPSILN*DTOR)*SIN(ALPHA*DTOR))/DTOR C Calculate Solar Hour Angle THSUN=(GHA360+ALON*RADDEG )*degrad SHA=GHA360+(ALON*RADDEG) c --- sun declination : SUNDEC = DEC*DEGRAD TRM111=SIN(ALAT)*SIN(DEC*DTOR) TRM112=-COS(ALAT)*COS(DEC*DTOR) TRM11=TRM111-TRM112 c --- solar noon altitude in degrees : SOLALT=ASIN(TRM11)/DTOR SUNBET = SOLALT c c --- cosine of the solar zenith angle : c coszen =sin(alat)*sin(sundec)+cos(alat)*cos(sundec)*cos(thsun) c if(coszen .le. 0.0) then coszen = 0.0 qatten = 0.0 else qatten = tau**(1./coszen) end if qzer = coszen * solar qdir = qzer * qatten qdiff = ((1.-aozone)*qzer - qdir) * 0.5 qtot = qdir + qdiff c c --- ( radiation as from Reed(1977), Simpson and Paulson(1979) ) c c c----------------------------------------------------------------------- c --- calculates the albedo as a function of the solar zenith angle : c --- ( after Payne jas 1972 ) c----------------------------------------------------------------------- c c --- solar zenith angle in degrees : c zen=raddeg*acos(coszen) c if(zen.ge.74.)then jab=.5*(90.-zen)+1. elseif(zen.ge.50.)then jab=.23*(74.-zen)+9. else jab=.10*(50.-zen)+15. endif c dzen=(za(jab)-zen)/dza(jab) c albedo=alb1(jab)+dzen*(alb1(jab+1)-alb1(jab)) c c --- calculates SHORT WAVE FLUX ( watt/m*m ) c --- ( Rosati,Miyakoda 1988 ; eq. 3.8 ) c c bb1=0.62 ! Reed bb2=0.636 ! Isemer et al. (1989) c qshort = qtot*(1. - bb1*acl + 0.0019*sunbet)*(1. - albedo) if (qshort.gt.qtot) qshort=qtot c 1000 return end