|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 !-------------------------------------- LICENCE BEGIN ------------------------------------ 00002 !Environment Canada - Atmospheric Science and Technology License/Disclaimer, 00003 ! version 3; Last Modified: May 7, 2008. 00004 !This is free but copyrighted software; you can use/redistribute/modify it under the terms 00005 !of the Environment Canada - Atmospheric Science and Technology License/Disclaimer 00006 !version 3 or (at your option) any later version that should be found at: 00007 !http://collaboration.cmc.ec.gc.ca/science/rpn.comm/license.html 00008 ! 00009 !This software is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; 00010 !without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 00011 !See the above mentioned License/Disclaimer for more details. 00012 !You should have received a copy of the License/Disclaimer along with this software; 00013 !if not, you can write to: EC-RPN COMM Group, 2121 TransCanada, suite 500, Dorval (Quebec), 00014 !CANADA, H9P 1J3; or send e-mail to service.rpn@ec.gc.ca 00015 !-------------------------------------- LICENCE END -------------------------------------- 00016 !copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%% 00017 ***S/P FLXSURF3 00018 * 00019 SUBROUTINE FLXSURF3BX(CMU, CTU, RIB, FTEMP, FVAP, ILMO, 00020 1 UE, FCOR, TA , QA , ZU, ZT, VA, 00021 1 TG , QG , H , Z0 , Z0T, 00022 1 LZZ0, LZZ0T, FM, FH, N ) 00023 ! 00024 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00025 USE PARKIND1 ,ONLY : JPRB 00026 ! 00027 #include "impnone.h" 00028 INTEGER N 00029 REAL CMU(N),CTU(N),RIB(N),FCOR(N),ILMO(N) 00030 REAL FTEMP(N),FVAP(N),TA(N),QA(N),ZU(N),ZT(N),VA(N) 00031 REAL TG(N),QG(N),H(N),Z0(N),UE(N) 00032 REAL Z0T(N),LZZ0(N),LZZ0T(N) 00033 REAL fm(N),fh(N) 00034 * 00035 *Author 00036 * Y.Delage (Jul 1990) 00037 *Revision 00038 * 001 G. Pellerin (Jun 94) New function for unstable case 00039 * 002 G. Pellerin (Jui 94) New formulation for stable case 00040 * 003 B. Bilodeau (Nov 95) Replace VK by KARMAN 00041 * 004 M. Desgagne (Dec 95) Add safety code in function ff 00042 * and ensures that RIB is non zero 00043 * 005 R. Sarrazin (Jan 96) Correction for H 00044 * 006 C. Girard (Nov 95) - Diffuse T instead of Tv 00045 * 007 G. Pellerin (Feb 96) Revised calculation for H (stable) 00046 * 008 G. Pellerin (Feb 96) Remove corrective terms to CTU 00047 * 009 Y. Delage and B. Bilodeau (Jul 97) - Cleanup 00048 * 010 Y. Delage (Feb 98) - Addition of HMIN 00049 * 011 D. Talbot and Y. Delage (Jan 02) - 00050 * Correct bug of zero divide by dg in loop 35 00051 * 012 Y. Delage (Oct 03) - Set top of surface layer at ZU +Z0 00052 * - Output UE instead of UE**2 and rename subroutine 00053 * - Change iteration scheme for stable case 00054 * - Introduce log-linear profile for near-neutral stable cases 00055 * - set VAMIN inside flxsurf and initialise ILMO and H 00056 * - Put stability functions into local functions via stabfunc.h 00057 * 013 Y. Delage (Sep 04) - Input of wind and temperature/humidity 00058 * at different levels 00059 * 014 R. McTaggart-Cowan and B. Bilodeau (May 2006) - 00060 * Clean up stabfunc.h 00061 * 015 L. Spacek (Dec 07) - Correction of the log-linear profile 00062 * Double precision for rib calculations 00063 * 00064 *Object 00065 * to calculate surface layer transfer coefficients and fluxes 00066 * 00067 *Arguments 00068 * 00069 * - Output - 00070 * CMU transfer coefficient of momentum times UE 00071 * CTU transfer coefficient of temperature times UE 00072 * RIB bulk Richardson number 00073 * FTEMP temperature flux 00074 * FVAP vapor flux 00075 * ILMO (1/length of Monin-Obukov) 00076 * UE friction velocity 00077 * H height of the boundary layer 00078 * FM momentum stability function 00079 * FH heat stability function 00080 * LZZ0 log ((zu+z0)/z0) 00081 * LZZ0T log ((zt+z0)/z0t) 00082 * 00083 * - Input - 00084 * FCOR Coriolis factor 00085 * ZU height of wind input (measured from model base at topo height + Z0) 00086 * ZT height of temperature and humidity input 00087 * TA potential temperature at ZT 00088 * QA specific humidity at ZT 00089 * VA wind speed at ZU 00090 * TG surface temperature 00091 * QG specific humidity at the surface 00092 * Z0 roughness length for momentum flux calculations 00093 * Z0T roughness length for heat/moisture flux calculations 00094 * N horizontal dimension 00095 * 00096 * 00097 #include "surfcon.h" 00098 #include "consphy.h" 00099 * 00100 ** 00101 * 00102 INTEGER J 00103 INTEGER IT,ITMAX 00104 REAL HMAX,CORMIN,EPSLN 00105 REAL CM,CT,ZP 00106 REAL F,G,DG 00107 REAL hi,HE,HS,unsl 00108 REAL*8 DTHV,TVA,TVS 00109 REAL HL,VAMIN,U 00110 REAL CS,XX,XX0,YY,YY0 00111 REAL ZB,DD,ILMOX 00112 REAL DF,ZZ,betsasx 00113 REAL aa,bb,cc 00114 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00115 SAVE HMAX,CORMIN,EPSLN 00116 SAVE ITMAX 00117 SAVE VAMIN 00118 * 00119 DATA CORMIN, HMAX /0.7E-4 , 1500.0/ 00120 DATA ITMAX / 3 / 00121 DATA VAMIN /0.1/ 00122 DATA EPSLN / 1.0e-05 / 00123 * 00124 DF(ZZ)=(1-ZZ*hi)*sqrt(1+4*AS*BETA*unsl*ZZ/(1-ZZ*hi)) 00125 CS=AS*2.5 00126 betsasx=1./asx 00127 * 00128 IF (LHOOK) CALL DR_HOOK('FLXSURF3BX',0,ZHOOK_HANDLE) 00129 * 00130 DO J=1,N 00131 LZZ0 (J)=1+ZU(J)/Z0(J) 00132 LZZ0T(J)=(ZT(J)+Z0(J))/Z0T(J) 00133 ENDDO 00134 * 00135 call vslog(LZZ0T,LZZ0T,N) 00136 call vslog(LZZ0 ,LZZ0 ,N) 00137 * 00138 DO J=1,N 00139 * 00140 * CALCULATE THE RICHARDSON NUMBER 00141 ZP=ZU(J)**2/(ZT(J)+Z0(J)-Z0T(J)) 00142 u=max(vamin,va(j)) 00143 tva=(1.d0+DELTA*QA(J))*TA(J) 00144 tvs=(1.d0+DELTA*QG(J))*TG(J) 00145 dthv=tva-tvs 00146 RIB(J)=GRAV/(tvs+0.5*dthv)*ZP*dthv/(u*u) 00147 if (rib(j).ge.0.0) rib(j) = max(rib(j), EPSLN) 00148 if (rib(j).lt.0.0) rib(j) = min(rib(j),-EPSLN) 00149 * 00150 * FIRST APPROXIMATION TO ILMO 00151 IF(RIB(J).GT.0.) THEN 00152 FM(J)=LZZ0(J)+CS*RIB(J)/max(2*z0(j),1.0) 00153 FH(J)=BETA*(LZZ0T(J)+CS*RIB(J))/ 00154 1 max(sqrt(z0(j)*z0t(j)),1.0) 00155 ILMO(J)=RIB(J)*FM(J)*FM(J)/(ZP*FH(J)) 00156 F=MAX(ABS(FCOR(J)),CORMIN) 00157 H(J)=BS*sqrt(KARMAN*u/(ILMO(J)*F*fm(j))) 00158 ELSE 00159 FM(J)=LZZ0(J)-min(0.7+log(1-rib(j)),LZZ0(J)-1) 00160 FH(J)=BETA*(LZZ0T(J)-min(0.7+log(1-rib(j)),LZZ0T(J)-1)) 00161 ILMO(J)=RIB(J)*FM(J)*FM(J)/(ZP*FH(J)) 00162 ENDIF 00163 ENDDO 00164 * 00165 * - - - - - - - - - BEGINNING OF ITERATION LOOP - - - - - - - - - - - 00166 DO 35 IT=1,ITMAX 00167 DO 35 J=1,N 00168 u=max(vamin,va(j)) 00169 ZP=ZU(J)**2/(ZT(J)+Z0(J)-Z0T(J)) 00170 IF(RIB(J).GT.0.) THEN 00171 *---------------------------------------------------------------------- 00172 * STABLE CASE 00173 ILMO(J)=max(EPSLN,ILMO(J)) 00174 hl=(ZU(J)+10*Z0(J))*FACTN 00175 F=MAX(ABS(FCOR(J)),CORMIN) 00176 hs=BS*sqrt(KARMAN*u/(ILMO(J)*F*fm(j))) 00177 H(J)=MAX(HMIN,hs,hl,factn/(4*AS*BETA*ilmo(j))) 00178 hi=1/h(j) 00179 *CDIR IEXPAND 00180 fm(J)=LZZ0(J)+psi(ZU(J)+Z0(J),hi,ilmo(j))-psi(Z0(J),hi,ilmo(j)) 00181 *CDIR IEXPAND 00182 fh(J)=BETA*(LZZ0T(J)+psi(ZT(J)+Z0(J),hi,ilmo(j))-psi(Z0T(J),hi, 00183 1 ilmo(j))) 00184 unsl=ILMO(J) 00185 DG=-ZP*FH(J)/(FM(J)*FM(J))*(1+beta*(DF(ZT(J)+Z0(J))-DF(Z0T(J)))/ 00186 1 (2*FH(J))-(DF(ZU(J)+Z0(J))-DF(Z0(J)))/FM(J)) 00187 *---------------------------------------------------------------------- 00188 * UNSTABLE CASE 00189 ELSE 00190 ILMO(J)=MIN(0.,ILMO(J)) 00191 *CDIR IEXPAND 00192 FM(J)=fmi(zu(j)+z0(j),z0 (j),lzz0 (j),ilmo(j),xx,xx0) 00193 *CDIR IEXPAND 00194 FH(J)=fhi(zt(j)+z0(j),z0t(j),lzz0t(j),ilmo(j),yy,yy0) 00195 DG=-ZP*FH(J)/(FM(J)*FM(J))*(1+beta/FH(J)*(1/YY-1/YY0)-2/FM(J)* 00196 1 (1/XX-1/XX0)) 00197 ENDIF 00198 *---------------------------------------------------------------------- 00199 IF(IT.LT.ITMAX) THEN 00200 G=RIB(J)-FH(J)/(FM(J)*FM(J))*ZP*ILMO(J) 00201 ILMO(J)=ILMO(J)-G/DG 00202 ENDIF 00203 35 CONTINUE 00204 * - - - - - - - - - END OF ITERATION LOOP - - - - - - - - - - - - - - 00205 * 00206 DO 80 J=1,N 00207 u=max(vamin,va(j)) 00208 if(asx.lt.as) then 00209 *---------------------------------------------------------------------- 00210 * CALCULATE ILMO AND STABILITY FUNCTIONS FROM LOG-LINEAR PROFILE 00211 * (SOLUTION OF A QUADRATIC EQATION) 00212 * 00213 zb=zu(j)/(zt(j)+z0(j)-z0t(j)) 00214 * DISCRIMINANT 00215 dd=(beta*lzz0t(j)*zb)**2-4*rib(j)*asx*lzz0(j)* 00216 1 (beta*lzz0t(j)*zb-lzz0(j)) 00217 if(rib(j).gt.0..and.rib(j).lt.betsasx.and.dd.ge.0.) then 00218 * COEFFICIENTS 00219 aa=asx*asx*rib(j)-asx 00220 bb=-beta*lzz0t(j)*zb+2*rib(j)*asx*lzz0(j) 00221 cc=rib(j)*lzz0(j)**2 00222 * SOLUTION 00223 if(bb>=0)then 00224 ilmox=(-bb-sqrt(dd)) 00225 1 /(2*zu(j)*aa) 00226 else 00227 ilmox=2*cc/(zu(j)*(-bb+sqrt(dd))) 00228 endif 00229 if(ilmox.lt.ilmo(j)) then 00230 ilmo(j)=ilmox 00231 fm(j)=lzz0(j)+asx*zu(j)*ilmox 00232 fh(j)=beta*lzz0t(j)+asx*(zt(j)+z0(j)-z0t(j))*ilmox 00233 endif 00234 endif 00235 endif 00236 *---------------------------------------------------------------------- 00237 CM=KARMAN/FM(J) 00238 CT=KARMAN/FH(J) 00239 UE(J)=u*CM 00240 CMU(J)=CM*UE(J) 00241 CTU(J)=CT*UE(J) 00242 if (rib(j).gt.0.0) then 00243 * stable case 00244 H(J)=MIN(H(J),hmax) 00245 else 00246 * unstable case 00247 F=MAX(ABS(FCOR(J)),CORMIN) 00248 he=max(HMIN,0.3*UE(J)/F) 00249 H(J)=MIN(he,hmax) 00250 endif 00251 FTEMP(J)=-CTU(J)*(TA(J)-TG(J)) 00252 FVAP(J)=-CTU(J)*(QA(J)-QG(J)) 00253 80 CONTINUE 00254 IF (LHOOK) CALL DR_HOOK('FLXSURF3BX',1,ZHOOK_HANDLE) 00255 CONTAINS 00256 #include "stabfunc2.h" 00257 END SUBROUTINE FLXSURF3BX
1.8.0
