C $Header: /u/gcmpack/MITgcm/pkg/cheapaml/cheapaml_lanl_flux.F,v 1.6 2012/12/28 23:25:50 jmc Exp $ C $Name: $ #include "CHEAPAML_OPTIONS.h" #undef ALLOW_THSICE CBOP C !ROUTINE: CHEAPAML_LANL_FLUX C !INTERFACE: SUBROUTINE CHEAPAML_LANL_FLUX( I i,j,bi,bj, O fsha, flha, evp, xolw, ssqt, q100 ) C !DESCRIPTION: C ================================================================== C SUBROUTINE cheapaml_LANL_flux C ================================================================== C o compute surface fluxes using LANL algorithms C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" c#include "FFIELDS.h" c#ifdef ALLOW_THSICE c#include "THSICE_VARS.h" c#endif #include "CHEAPAML.h" C !INPUT PARAMETERS: INTEGER i,j,bi,bj C !OUTPUT PARAMETERS: _RL fsha, flha, evp, xolw, ssqt, q100 CEOP C Output: C ustress, vstress - wind stress C fsha - sensible heat flux C flha - latent heat flux C xolw - oceanic upwelled long wave radiation C ssqt - sat. specific humidity at atm layer top C Input: C uwind, vwind - mean wind speed (m/s) C Tair - mean air temperature (K) at height ht (m) C theta(k=1) - sea surface temperature (C) C Qair - Specific humidity kg/kg C Solar - short wave net solar flux at surface (W/m^2) C Tr - relaxation profile for temperature on boundaries (C) C qr - relaxation profile for specific humidity (kg/kg) C i,j,bi,bj - indices of data C !LOCAL VARIABLES: C iceornot :: variables to include seaice effect INTEGER iceornot _RL deltaTm _RL uss,usm,uw,vw _RL cheapaml_BulkCdn _RL to _RL t _RL t0,QaR _RL ssq, q _RL deltap, delq, pt, psx100, z100ol _RL rdn,ren,rhn,zice,zref _RL rd,re,rh,tta,ttas,toa,ttt _RL ustar,tstar,qstar,ht,hu,hq _RL aln,cdalton,czol,psim_fac _RL huol,stable,xsq,x,psimh,psixh _RL clha, csha INTEGER niter_bulk,iter C useful values C hardwire atmospheric relative humidity at 80% QaR=0.8 _d 0 C factor to compute rainfall from specific humidity C inverse of time step deltaTm=1. _d 0/deltaT C reference values to compute turbulent flux ht=zt hq=zq hu=zu zref = zt zice=.0005 _d 0 aln = log(ht/zref) C for iterating on turbulence niter_bulk = 5 cdalton = 0.0346000 _d 0 czol = zref*xkar*gravity psim_fac=5. _d 0 C determine wind stress IF(.NOT.useStressOption)THEN if (maskC(i,j,1,bi,bj).ne.0. _d 0) then #ifdef ALLOW_THSICE if (ICEMASK(i,j,bi,bj).gt.0. _d 0) then if (snowheight(i,j,bi,bj).gt.3. _d -1) then iceornot=2 else iceornot=1 endif else iceornot=0 endif #else iceornot=0 #endif uw=uwind(i,j,bi,bj) vw=vwind(i,j,bi,bj) uss=sqrt(uw**2+vw**2) usm=max(uss,1. _d 0) cheapaml_BulkCdn = cdrag_1/usm + cdrag_2 + cdrag_3*usm ustress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*uw vstress(i,j,bi,bj)= rhoa*cheapaml_BulkCdn*uss*vw else usm=0. _d 0 ustress(i,j,bi,bj) = 0. _d 0 vstress(i,j,bi,bj) = 0. _d 0 endif C wind stress computed ENDIF C diabatic and freshwater flux forcing to=theta(i,j,1,bi,bj) t=Tair(i,j,bi,bj) toa=to+Celsius2K tta=t+Celsius2K ttas=tta+gamma_blk*zref ttt=tta-( cheaphgrid(i,j,bi,bj)- zref)*gamma_blk pt=p0*(1-gamma_blk*cheaphgrid(i,j,bi,bj)/ttas) & **(gravity/gamma_blk/gasR) C specific humidities ssq= ssq0*exp( lath*(ssq1-ssq2/toa) ) / p0 ssqt = ssq0*exp( lath*(ssq1-ssq2/ttt) ) / pt C saturation no more at the top: ssqt = 0.7 _d 0*ssq if (useFreshWaterFlux) then q=qair(i,j,bi,bj) else q=QaR * ssq endif C adjust temperature from reference height to formula height deltap = t - to + gamma_blk*(zref-ht) delq = q - ssq ttas = tta+gamma_blk*(zref-ht) t0 = ttas*(1. _d 0 + humid_fac*q) C initialize estimate exchange coefficients rdn=xkar/(log(zref/zice)) rhn=rdn ren=rdn C calculate turbulent scales ustar=rdn*usm tstar=rhn*deltap qstar=ren*delq C iteration with psi-functions to find transfer coefficients do iter=1,niter_bulk huol = czol/ustar**2 *(tstar/t0 + & qstar/(1. _d 0/humid_fac+q)) huol = sign( min(abs(huol),10. _d 0), huol) stable = 5. _d -1 + sign(5. _d -1 , huol) xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*huol)),1. _d 0) x = sqrt(xsq) psimh = -5. _d 0*huol*stable + (1. _d 0-stable)* & (2. _d 0*log(5. _d -1*(1. _d 0+x)) + & 2. _d 0*log(5. _d -1*(1. _d 0+xsq)) - & 2. _d 0*atan(x) + pi*.5 _d 0) psixh = -5. _d 0*huol*stable + (1. _d 0-stable)* & (2. _d 0*log(5. _d -1*(1. _d 0+xsq))) C Update the transfer coefficients rd = rdn/(1. _d 0 + rdn*(aln-psimh)/xkar) rh = rhn/(1. _d 0 + rhn*(aln-psixh)/xkar) re = rh C Update ustar, tstar, qstar using updated, shifted coefficients. ustar = rd*usm qstar = re*delq tstar = rh*deltap enddo usm=max(uss,0.5 _d 0) csha = rhoa*cpair*usm*rh*rd clha = rhoa*lath*usm*re*rd fsha = csha*deltap flha = clha*delq evp = -flha/lath C the sensible and latent heat fluxes, fsha and flha, C are computed so that positive values are downward. C the convention for cheapaml is upward fluxes are positive, C so they must be multiplied by -1 fsha=-fsha flha=-flha C oceanic upwelled long wave xolw=stefan*(toa)**4 C compute specific humidity at 100m huol = czol/ustar**2 *(tstar/t0 + & qstar/(1. _d 0/humid_fac+q)) huol = sign( min(abs(huol),10. _d 0), huol) stable = 5. _d -1 + sign(5. _d -1 , huol) z100ol = 100. _d 0 *xkar*gravity/ustar**2 *(tstar/t0 & + qstar/(1. _d 0/humid_fac+q)) xsq = max(sqrt(abs(1. _d 0 - 16. _d 0*z100ol)),1. _d 0) x = sqrt(xsq) psx100 = -5. _d 0*z100ol*stable + (1. _d 0-stable)* & (2. _d 0*log(5. _d -1*(1. _d 0+xsq))) q100=ssq+qstar*(dlog(100. _d 0/zice)-psx100) RETURN END