C A DEMONSTRATION PROGRAM FOR CALCULATION OF ATMOSPHERIC C SURFACE LAYER FLUXES IN CONVECTIVE CONDITIONS C USING TURBULENCE CLOSURE MODELS C C AUTHOR: LECH LOBOCKI, Warsaw University of Technology C Nowowiejska 20, 00-653 Warsaw, Poland C E: llobocki@is.pw.edu.pl C C SCIENTIFIC DOCUMENTATION: C Lech Lobocki, 2001: "Calculation of surface fluxes under C convective conditions by turbulence closure models", C Journal of Applied Meteorology, 40: $-$ (JAM 1690). C Lech Lobocki, 2001: "An explicit agorithm for calculating C surface-layer parameters in convective conditions derived C from a turbulence closure model", C Journal of Applied Meteorology, 40: $-$ (JAM 1697). C C VERSION: 1.00, December 30, 2000 C C LANGUAGE: FORTRAN-77 C C TERMS OF USE: GPL LICENSE (see below). C In addition, the use of this code must be acknowledged C in scientific publications, project reports, or C any other documents describing research or other work C where this code was used, either directly or as a coding example. C C BUGS: C So far, none reported. C Please direct bug reports and comments via e-mail to the address C given above. C C NOTES: C C 1) This code is intended mainly as an illustration of methods C presented in scientific publications. Therefore, it aims C at clarity and convertibility rather than at computational C efficiency. For massive computations, it must be rewritten C into a form suitable for a given machine and master-code C architectures. C 2) This file contains four sections: section 1 is a main program C that tabulates MF(dU,dT) and HF(dU,dT) functions; section 2 C implements a particular turbulence closure model (here, the C Mellor-Yamada Level 2); section 3 is an iterative 'approximate' C algorithm described in JAM1690; and section 4 is an explicit C algorithm presented in JAM1697. C 3) Inline comments refer to equations numbers in JAM1690 (section 3) C and JAM1697 (section 4) C C GPL LICENSING INFO: C Copyright (C) 2001 Lech Lobocki C C This program is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the License, or C (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public License C along with this program; if not, write to the Free Software C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. C C ---------------------------------------------------------------- C C SECTION 1: Main program that produces arrays of fluxes C used to prepare Fig. 9 in JAM1690 C C ---------------------------------------------------------------- C parameter (N=60, L=50) double precision MF,HF,dU,dT,U,ddT,z0,zSL,ro,cp,uw,wT double precision AMINARR,AMAXARR double precision dTMin, dTMax, Umin, Umax dimension MF(N,L),HF(N,L) ro=1.27D0 cp=1.005D3 zSL=1D1 z0=1D-2 Umin=0D0 Umax=5D0 dTmin=-2D1 dTmax=0D0 dU=(Umax-Umin)/DFLOAT(L-1) ddT=(dTmax-dTmin)/DFLOAT(N-1) do i=1,N dT=dTmin+(i-1)*ddT do j=1,L U=Umin+(j-1)*dU call SLFX (z0,zSL,U,dT,wT,uw) C or, use SLFXA, if you prefer C call SLFXA (z0,zSL,U,dT,wT,uw) MF(i,j)=uw*ro HF(i,j)=wT*ro*cp enddo enddo OPEN (1,FILE='stress.grd') CALL GRDWRITE (1, n, l, dTmin, dTmax, Umin, Umax, & AMINARR(MF,N*L), AMAXARR(MF,N*L), MF) CLOSE (1) OPEN (1,FILE='hflux.grd') CALL GRDWRITE (1, n, l, dTmin, dTmax, Umin, Umax, & AMINARR(HF,N*L), AMAXARR(HF,N*L), HF) CLOSE (1) stop end SUBROUTINE GRDWRITE (LUN, IM, JM, XMIN, XMAX, & YMIN, YMAX, ZMIN, ZMAX, A) double precision A(IM,JM),xmin,xmax,ymin,ymax,zmin,zmax WRITE (LUN, '(A4)') 'DSAA' WRITE (LUN, 301) IM,JM WRITE (LUN, 303) XMIN, XMAX WRITE (LUN, 303) YMIN, YMAX WRITE (LUN, 303) ZMIN, ZMAX DO J=1,JM WRITE (LUN, '(10(G12.5,1X))') (A(I,J),I=1,IM) ENDDO RETURN 301 FORMAT (2(I4,1X)) 303 FORMAT (2(G12.5,1X)) END double precision FUNCTION AMINARR (A,N) double precision A(N),F F=A(1) DO I=2,N F=DMIN1(F,A(I)) ENDDO AMINARR=F RETURN END double precision FUNCTION AMAXARR (A,N) double precision A(N),F F=A(1) DO I=2,N F=DMAX1(F,A(I)) ENDDO AMAXARR=F RETURN END C ---------------------------------------------------------------- C C SECTION 2: Turbulence closure model imlementation. C Here, the Mellor-Yamada Level 2 model is used. C For discussion of values of the model constants, C see JAM1690. C You may replace this section by a model of your own... C C ---------------------------------------------------------------- block data double precision vKarm, BuyPar double precision A1,A2,B1,B2,C1 common /physcn/ vKarm, BuyPar common /MY2/ A1,A2,B1,B2,C1 data vKarm, BuyPar /0.4D0, 0.034D0/ data A1,A2,B1,B2,C1 /0.69D0,0.57D0,16.6D0,7.9D0, 0.061D0/ end double precision function M_fhneu() M_fhneu=0.92D0 return end double precision function M_psi_m_apx(dzeta) double precision dzeta if (dzeta.LT.-1.0D0) stop 'BAD DZETA' M_psi_m_apx=((0.3889D0*dzeta+1.2488D0)*dzeta+2.0358D0)*dzeta ! Eq. (83) return end double precision function M_psi_h_apx(dzeta) double precision dzeta if (dzeta.LT.-1.0D0) stop 'BAD DZETA' M_psi_h_apx=((0.6783D0*dzeta+1.8673D0)*dzeta+2.5408D0)*dzeta ! Eq. (84) return end double precision function M_SH_inf() double precision A1,A2,B1,B2,C1 common /MY2/ A1,A2,B1,B2,C1 M_SH_inf=DSQRT(2.0D0)*A2/B1*(12D0*A1+B1+3D0*B2) ! Table 2 return end double precision function M_SM_inf() double precision A1,A2,B1,B2,C1 double precision M_SH_inf common /MY2/ A1,A2,B1,B2,C1 M_SM_inf=M_SH_inf()*A1/A2*(B1-3D0*B1*C1+12D0*A1+9D0*A2) & /(B1+3D0*(B2+A1)) ! Table 2 return end double precision function M_apx1(x,z_z0) double precision x,z_z0 c 1st approximation step after normalization by dividing c through ln(z/z0)*(1-RiB)**-0.25 double precision aaa,bbb,Om_inf,B1_1_3,onethird double precision A1,A2,B1,B2,C1 double precision M_SH_inf double precision M_SM_inf common /MY2/ A1,A2,B1,B2,C1 onethird=1.0D0/3.0D0 B1_1_3=B1**onethird aaa=A2*B1_1_3*(1.0D0-6D0*A1/B1)*DLOG(z_z0) aaa=aaa**(-4) ! Eq. (19), JAM1697 Om_inf=3D0*DSQRT(2.0D0)*M_SH_inf()/((M_SM_inf())**2*B1_1_3) & *((z_z0**onethird)-1D0) bbb=Om_inf**(-3) ! Eq. (20), JAM1697 M_apx1=(aaa-bbb*x)**(-0.25) ! Eq. (18), JAM1697 return end double precision function M_OmegaN_apx(RiB,z_z0) double precision RiB,z_z0 double precision lambda, alpha, AA, x, x1, x0, residuum double precision M_apx1 lambda=DLOG10(z_z0) AA= 0.015D0 * lambda + 0.05D0 ! Eq. (23), JAM1697 alpha=-0.125D0 * lambda + 1.25D0 ! Eq. (24), JAM1697 x1=0.55D0 * lambda -1.6D0 ! Eq. (25), JAM1697 x0= 1.125D0 * lambda -3.25D0 ! Eq. (26), JAM1697 x=DLOG10(-RiB) residuum=AA*(x-x0)*DEXP(-(alpha*(x-x1))**2); ! Eq. (22), JAM1697 M_OmegaN_apx=M_apx1(RiB,z_z0)/DLOG(z_z0)+residuum ! Eq. (21), JAM1697 return end double precision function M_CFlux(dT, z0, zSL) ! returns free-convection value of the kinematic heat flux given ! the temperature difference dT between two levels, z0 and zSL double precision dz13, onethird double precision dT, z0, zSL double precision A1,A2,B1,B2,C1 double precision vKarm, BuyPar common /physcn/ vKarm, BuyPar common /MY2/ A1,A2,B1,B2,C1 if (dT.GE.0.0D0) STOP 'BAD DT' onethird=1.0D0/3.0D0 dz13=z0**(-onethird)-zSL**(-onethird) M_CFlux=vKarm*vKarm*DSQRT(BuyPar)/ & B1*((-dT*A2*(B1+12.0D0*A1+3.0D0*B2)/(3.0D0*dz13))**1.5D0) ! Eq. (87) return end double precision function M_Imc_apx(zSL,mL_z) double precision zSL,mL_z, onethird double precision M_IMc_c_inf onethird=1.0D0/3.0D0 if (mL_z.EQ.0D0) then M_Imc_apx=M_IMc_c_inf()*(zSL**(-onethird)) else M_Imc_apx=M_IMc_c_inf()*(zSL**(-onethird)) & *(1D0-(mL_z**onethird))*(mL_z**(-onethird)) ! Eq. (92) endif return end double precision function M_Ihc_apx(zSL,mL_z) double precision zSL,mL_z, onethird double precision M_IHc_c_inf onethird=1.0D0/3.0D0 if (mL_z.EQ.0D0) then M_Ihc_apx=M_IHc_c_inf()*(zSL**(-onethird)) else M_Ihc_apx=M_IHc_c_inf()*(zSL**(-onethird)) & *(1D0-(mL_z**onethird))*(mL_z**(-onethird)) ! Eq. (92) endif return end double precision function M_Imc_c_inf() double precision vKarm, BuyPar double precision A1,A2,B1,B2,C1 double precision M_IHc_c_inf common /physcn/ vKarm, BuyPar common /MY2/ A1,A2,B1,B2,C1 M_Imc_c_inf=M_IHc_c_inf()*A2/A1*(B1+3D0*(B2+A1)) & /(B1-3D0*B1*C1+12D0*A1+9D0*A2) ! Eq. (95) return end double precision function M_Ihc_c_inf() double precision vKarm, BuyPar double precision A1,A2,B1,B2,C1 common /physcn/ vKarm, BuyPar common /MY2/ A1,A2,B1,B2,C1 M_Ihc_c_inf=3.0D0*(B1/(vKarm*vKarm))**(2.0D0/3.0D0) & /(A2*(B1+12D0*A1+3D0*B2)) return ! Eq. (93) end C ---------------------------------------------------------------- C C SECTION 3: iterative 'approximate' algorithm described in JAM1690. C C ---------------------------------------------------------------- subroutine SLFX (z0,zSL,U,dT,wT,uw) C C ---------------------------------------------------------------- C C INPUT PARAMETERS: C z0: lower computational level (not necessarily aerodynamic roughness) C zSL: upper computational level (assumed somewhere in the ASL) C U: wind speed difference between these two levels C dT: potential temperature difference between these two levels C OUTPUT PARAMETERS: C uw: kinematic vertical turbulent momentum flux C wT: kinematic vertical turbulent sensible heat flux C UNITS: SI C C ---------------------------------------------------------------- C double precision z0,zSL,U,dT,wT,uw double precision L_inv_old, mL_z, mL_z0, dzeta, dzeta0 double precision px, pL, L, L_inv, onethird, lnz_z0, lnL_z0 double precision ustar, Tstar double precision IMd, IMc, IHd, IHc double precision eps double precision M_CFlux, M_psi_m_apx, M_psi_h_apx double precision M_fhneu double precision vKarm, BuyPar double precision M_Ihc_apx, M_Imc_apx common /physcn/ vKarm, BuyPar data eps /1.0D-5/ if (U.LE.0) then if (dT.EQ.0) then ! no shear, no stratification => no turbulence wT=0.0 uw=0.0 else ! free convection uw=0.0 wT=M_CFlux(dT,z0,zSL) ! Eq. (87) endif else if (dT.EQ.0.0D0) then ! neutral wT=0 uw=(vKarm*U/DLOG(zSL/z0))**2 else onethird=1D0/3D0 lnz_z0=DLOG(zSL/z0) ! take a neutral asymptotic value for uw as the first guess uw=(vKarm*U/DLOG(zSL/z0))**2 ! take convective asymptotic value for wT as the first guess wT=M_CFlux(dT, z0, zSL) ! Eq. (87) ! and use the Monin-Obukhov lenght scale resulting from those ! as the first guess L_inv=-(vKarm*BuyPar*wT)/(uw**1.5D0) ! Eq. (3) L_inv_old=L_inv*(1D0+2D0*eps) ! just to start DO WHILE ((DABS(L_inv_old-L_inv)).GT.(eps*DABS(L_inv))) L_inv_old=L_inv L=1/L_inv if (z0.GT.DABS(L)) then ! near free convection, use convective scaling only mL_z=-L/zSL mL_z0=-L/z0 px=(BuyPar*wT)**(-onethird) wT=-dT/(px*(M_IHc_apx(zSL,mL_z)-M_IHc_apx(z0,mL_z0))) ! Eq. (97) uw=U/(px*(M_IMc_apx(zSL,mL_z)-M_IMc_apx(z0,mL_z0))) ! Eq. (97) else if (zSL.LT.DABS(L)) then ! near-neutral, use dynamic scaling dzeta=zSL/L dzeta0=z0/L ustar=U*vKarm/(lnz_z0+M_psi_M_apx(dzeta) & -M_psi_M_apx(dzeta0)) ! Eq. (41) Tstar=dT*vKarm/(M_fhneu()*(lnz_z0+M_psi_H_apx(dzeta) & -M_psi_H_apx(dzeta0))) ! Eq. (42) uw=ustar*ustar ! Eq. (1) wT=-ustar*Tstar ! Eq. (2) else ! z0<|L| no turbulence wT=0.0 uw=0.0 else ! free convection uw=0.0 wT=M_CFlux(dT,z0,zSL) endif else if (dT.EQ.0.0D0) then ! neutral wT=0 uw=(vKarm*U/DLOG(zSL/z0))**2 else onethird=1D0/3D0 lnz_z0=DLOG(zSL/z0) RiB_=(BuyPar*(zSL-z0)*dT)/(U*U) L=(zSL-z0)/(RiB_*M_OmegaN_apx(RiB_,zSL/z0)*lnz_z0) if (z0.GT.DABS(L)) then ! near free convection, use convective scaling only } mL_z=-L/zSL mL_z0=-L/z0 wT=DSQRT(BuyPar)*((-dT/(M_IHc_apx(zSL,mL_z) & -M_IHc_apx(z0,mL_z0)))**1.5) px=(BuyPar*wT)**(-onethird) uw=U/(px*(M_IMc_apx(zSL,mL_z)-M_IMc_apx(z0,mL_z0))) else if (zSL.LT.DABS(L)) then ! near-neutral, use dynamic scaling dzeta=zSL/L dzeta0=z0/L ustar=U*vKarm/(lnz_z0+M_psi_M_apx(dzeta) & -M_psi_M_apx(dzeta0)) Tstar=dT*vKarm/(M_fhneu()*(lnz_z0+M_psi_H_apx(dzeta) & -M_psi_H_apx(dzeta0))) uw=ustar*ustar wT=-ustar*Tstar else ! z0<|L|