[1] | 1 | !---------------------------------------------------------------------- |
---|
| 2 | ! |
---|
| 3 | MODULE MODULE_SF_QNSESFC |
---|
| 4 | ! |
---|
| 5 | !---------------------------------------------------------------------- |
---|
| 6 | ! |
---|
| 7 | USE MODULE_MODEL_CONSTANTS |
---|
| 8 | USE MODULE_DM, ONLY : WRF_DM_MAXVAL |
---|
| 9 | ! |
---|
| 10 | !---------------------------------------------------------------------- |
---|
| 11 | ! |
---|
| 12 | ! REFERENCES: Janjic (2002), NCEP Office Note 437 |
---|
| 13 | ! |
---|
| 14 | ! ABSTRACT: |
---|
| 15 | ! MYJSFC GENERATES THE SURFACE EXCHANGE COEFFICIENTS FOR VERTICAL |
---|
| 16 | ! TURBULENT EXCHANGE BASED UPON MONIN_OBUKHOV THEORY WITH |
---|
| 17 | ! VARIOUS REFINEMENTS. |
---|
| 18 | ! |
---|
| 19 | !---------------------------------------------------------------------- |
---|
| 20 | ! |
---|
| 21 | INTEGER :: ITRMX=5 ! Iteration count for sfc layer computations |
---|
| 22 | ! |
---|
| 23 | REAL,PARAMETER :: EPSQ2L=0.01 |
---|
| 24 | REAL,PARAMETER :: VKARMAN=0.4 |
---|
| 25 | REAL,PARAMETER :: CAPA=R_D/CP,ELOCP=2.72E6/CP,RCAP=1./CAPA |
---|
| 26 | REAL,PARAMETER :: GOCP02=G/CP*2.,GOCP10=G/CP*10. |
---|
| 27 | ! REAL,PARAMETER :: EPSU2=1.E-4,EPSUST=0.07,EPSZT=1.E-5 |
---|
| 28 | ! ECMWF sets lower limit on ln(Z0h)=-20 --> EPSZT=2.E-9 |
---|
| 29 | REAL,PARAMETER :: EPSU2=1.E-6,EPSUST=1.E-9,EPSZT=1.E-28 |
---|
| 30 | REAL,PARAMETER :: A2S=17.2693882,A3S=273.16,A4S=35.86 |
---|
| 31 | REAL,PARAMETER :: SEAFC=0.98,PQ0SEA=PQ0*SEAFC |
---|
| 32 | REAL,PARAMETER :: BETA=1./273.,CZIL=0.1,EXCML=0.0001,EXCMS=0.0001 & |
---|
| 33 | !old REAL,PARAMETER :: BETA=1./273.,CZIL=0.1,EXCML=0.001,EXCMS=0.001 & |
---|
| 34 | & ,GLKBR=10.,GLKBS=30.,PI=3.1415926 & |
---|
| 35 | & ,QVISC=2.1E-5,RIC=0.505,SMALL=0.35 & |
---|
| 36 | & ,SQPR=0.84,SQSC=0.84,SQVISC=258.2,TVISC=2.1E-5 & |
---|
| 37 | & ,USTC=0.7,USTR=0.225,VISC=1.5E-5,FH=1.01 & |
---|
| 38 | & ,WWST=1.2,ZTFC=0.1,TOPOFAC=9.0E-6 |
---|
| 39 | ! |
---|
| 40 | REAL,PARAMETER :: BTG=BETA*G,CZIV=SMALL*GLKBS & |
---|
| 41 | & ,GRRS=GLKBR/GLKBS & |
---|
| 42 | & ,RTVISC=1./TVISC,RVISC=1./VISC & |
---|
| 43 | & ,ZQRZT=SQSC/SQPR & |
---|
| 44 | & ,FZQ1=RTVISC*QVISC*ZQRZT & |
---|
| 45 | & ,FZQ2=RTVISC*QVISC*ZQRZT & |
---|
| 46 | & ,FZT1=RVISC *TVISC*SQPR & |
---|
| 47 | & ,FZT2=CZIV*GRRS*TVISC*SQPR & |
---|
| 48 | & ,FZU1=CZIV*VISC & |
---|
| 49 | & ,PIHF=0.5*PI & |
---|
| 50 | & ,PRT0=0.72 & |
---|
| 51 | & ,RQVISC=1./QVISC & |
---|
| 52 | & ,USTFC=0.018/G & |
---|
| 53 | & ,WWST2=WWST*WWST & |
---|
| 54 | & ,ZILFC=-CZIL*VKARMAN*SQVISC |
---|
| 55 | ! |
---|
| 56 | !---------------------------------------------------------------------- |
---|
| 57 | INTEGER, PARAMETER :: KZTM=10001,KZTM2=KZTM-2 |
---|
| 58 | ! |
---|
| 59 | REAL :: DZETA1,DZETA2,FH01,FH02,ZTMAX1,ZTMAX2,ZTMIN1,ZTMIN2 |
---|
| 60 | ! |
---|
| 61 | REAL,DIMENSION(KZTM) :: PSIH1,PSIH2,PSIM1,PSIM2 |
---|
| 62 | ! |
---|
| 63 | !---------------------------------------------------------------------- |
---|
| 64 | ! |
---|
| 65 | CONTAINS |
---|
| 66 | ! |
---|
| 67 | !---------------------------------------------------------------------- |
---|
| 68 | SUBROUTINE QNSESFC(ITIMESTEP,HT,DZ & |
---|
| 69 | & ,PMID,PINT,TH,T,QV,QC,U,V,Q2 & |
---|
| 70 | & ,TSK,QSFC,THZ0,QZ0,UZ0,VZ0 & |
---|
| 71 | & ,LOWLYR,XLAND & |
---|
| 72 | & ,USTAR,ZNT,Z0BASE,PBLH,MAVAIL,RMOL & |
---|
| 73 | & ,AKHS,AKMS & |
---|
| 74 | & ,RIB & |
---|
| 75 | & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC & |
---|
| 76 | & ,QGH,CPM,CT & |
---|
| 77 | & ,U10,V10,TSHLTR,TH10,QSHLTR,Q10,PSHLTR & |
---|
| 78 | & ,IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 79 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 80 | & ,ITS,ITE,JTS,JTE,KTS,KTE) |
---|
| 81 | !---------------------------------------------------------------------- |
---|
| 82 | ! |
---|
| 83 | IMPLICIT NONE |
---|
| 84 | ! |
---|
| 85 | !---------------------------------------------------------------------- |
---|
| 86 | INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 87 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 88 | & ,ITS,ITE,JTS,JTE,KTS,KTE |
---|
| 89 | ! |
---|
| 90 | INTEGER,INTENT(IN) :: ITIMESTEP |
---|
| 91 | ! |
---|
| 92 | INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: LOWLYR |
---|
| 93 | ! |
---|
| 94 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: HT,MAVAIL,TSK & |
---|
| 95 | & ,XLAND,Z0BASE |
---|
| 96 | ! |
---|
| 97 | REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: DZ & |
---|
| 98 | & ,PMID,PINT & |
---|
| 99 | & ,Q2,QC,QV & |
---|
| 100 | & ,T,TH & |
---|
| 101 | & ,U,V |
---|
| 102 | ! |
---|
| 103 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: FLX_LH,HFX,PSHLTR & |
---|
| 104 | & ,QFX,Q10,QSHLTR & |
---|
| 105 | & ,TH10,TSHLTR & |
---|
| 106 | & ,U10,V10 |
---|
| 107 | ! |
---|
| 108 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: AKHS,AKMS & |
---|
| 109 | & ,PBLH,QSFC |
---|
| 110 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: RIB |
---|
| 111 | ! |
---|
| 112 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: QZ0,RMOL,THZ0 & |
---|
| 113 | & ,USTAR,UZ0,VZ0 & |
---|
| 114 | & ,ZNT |
---|
| 115 | ! |
---|
| 116 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CHS,CHS2,CQS2 & |
---|
| 117 | & ,CPM,CT,FLHC,FLQC & |
---|
| 118 | & ,QGH |
---|
| 119 | !---------------------------------------------------------------------- |
---|
| 120 | !*** |
---|
| 121 | !*** LOCAL VARIABLES |
---|
| 122 | !*** |
---|
| 123 | INTEGER :: I,J,K,KFLIP,LMH,LPBL,NTSD |
---|
| 124 | ! |
---|
| 125 | REAL :: A,APESFC,B,BTGX,CWMLOW & |
---|
| 126 | & ,DQDT,DTDIF,DTDT,DUDT,DVDT & |
---|
| 127 | & ,FIS & |
---|
| 128 | & ,P02P,P10P,PLOW,PSFC,PTOP,QLOW,QS02,QS10 & |
---|
| 129 | & ,RAPA,RAPA02,RAPA10,RATIOMX,RDZ,SEAMASK,SM & |
---|
| 130 | & ,T02P,T10P,TEM,TH02P,TH10P,THLOW,THELOW,THM & |
---|
| 131 | & ,TLOW,TZ0,ULOW,VLOW,ZSL |
---|
| 132 | ! |
---|
| 133 | REAL,DIMENSION(KTS:KTE) :: CWMK,PK,Q2K,QK,THEK,THK,TK,UK,VK |
---|
| 134 | ! |
---|
| 135 | REAL,DIMENSION(KTS:KTE-1) :: EL,ELM |
---|
| 136 | ! |
---|
| 137 | REAL,DIMENSION(KTS:KTE+1) :: ZHK |
---|
| 138 | ! |
---|
| 139 | REAL,DIMENSION(ITS:ITE,JTS:JTE) :: THSK |
---|
| 140 | ! |
---|
| 141 | REAL,DIMENSION(ITS:ITE,KTS:KTE+1,JTS:JTE) :: ZINT |
---|
| 142 | ! |
---|
| 143 | !---------------------------------------------------------------------- |
---|
| 144 | !********************************************************************** |
---|
| 145 | !---------------------------------------------------------------------- |
---|
| 146 | ! |
---|
| 147 | !*** MAKE PREPARATIONS |
---|
| 148 | ! |
---|
| 149 | !---------------------------------------------------------------------- |
---|
| 150 | DO J=JTS,JTE |
---|
| 151 | DO K=KTS,KTE+1 |
---|
| 152 | DO I=ITS,ITE |
---|
| 153 | ZINT(I,K,J)=0. |
---|
| 154 | ENDDO |
---|
| 155 | ENDDO |
---|
| 156 | ENDDO |
---|
| 157 | ! |
---|
| 158 | DO J=JTS,JTE |
---|
| 159 | DO I=ITS,ITE |
---|
| 160 | ZINT(I,KTE+1,J)=HT(I,J) ! Z at bottom of lowest sigma layer |
---|
| 161 | PBLH(I,J)=-1. |
---|
| 162 | ! |
---|
| 163 | !!!!!!!!! |
---|
| 164 | !!!!!! UNCOMMENT THESE LINES IF USING ETA COORDINATES |
---|
| 165 | !!!!!!!!! |
---|
| 166 | !!!!!! ZINT(I,KTE+1,J)=1.E-4 ! Z of bottom of lowest eta layer |
---|
| 167 | !!!!!! ZHK(KTE+1)=1.E-4 ! Z of bottom of lowest eta layer |
---|
| 168 | ! |
---|
| 169 | ENDDO |
---|
| 170 | ENDDO |
---|
| 171 | ! |
---|
| 172 | DO J=JTS,JTE |
---|
| 173 | DO K=KTE,KTS,-1 |
---|
| 174 | KFLIP=KTE+1-K |
---|
| 175 | DO I=ITS,ITE |
---|
| 176 | ZINT(I,K,J)=ZINT(I,K+1,J)+DZ(I,KFLIP,J) |
---|
| 177 | ENDDO |
---|
| 178 | ENDDO |
---|
| 179 | ENDDO |
---|
| 180 | ! |
---|
| 181 | NTSD=ITIMESTEP |
---|
| 182 | ! |
---|
| 183 | IF(NTSD==1)THEN |
---|
| 184 | DO J=JTS,JTE |
---|
| 185 | DO I=ITS,ITE |
---|
| 186 | USTAR(I,J)=0.1 |
---|
| 187 | FIS=HT(I,J)*G |
---|
| 188 | SM=XLAND(I,J)-1. |
---|
| 189 | !!! Z0(I,J)=SM*Z0SEA+(1.-SM)*(Z0(I,J)*Z0MAX+FIS*FCM+Z0LAND) |
---|
| 190 | !!! ZNT(I,J)=SM*Z0SEA+(1.-SM)*(ZNT(I,J)*Z0MAX+FIS*FCM+Z0LAND) |
---|
| 191 | ENDDO |
---|
| 192 | ENDDO |
---|
| 193 | ENDIF |
---|
| 194 | ! |
---|
| 195 | !!!! IF(NTSD==1)THEN |
---|
| 196 | DO J=JTS,JTE |
---|
| 197 | DO I=ITS,ITE |
---|
| 198 | CT(I,J)=0. |
---|
| 199 | ENDDO |
---|
| 200 | ENDDO |
---|
| 201 | !!!! ENDIF |
---|
| 202 | ! |
---|
| 203 | !---------------------------------------------------------------------- |
---|
| 204 | setup_integration: DO J=JTS,JTE |
---|
| 205 | !---------------------------------------------------------------------- |
---|
| 206 | ! |
---|
| 207 | DO I=ITS,ITE |
---|
| 208 | ! |
---|
| 209 | !*** LOWEST LAYER ABOVE GROUND MUST BE FLIPPED |
---|
| 210 | ! |
---|
| 211 | LMH=KTE-LOWLYR(I,J)+1 |
---|
| 212 | ! |
---|
| 213 | PTOP=PINT(I,KTE+1,J) ! KTE+1=KME |
---|
| 214 | PSFC=PINT(I,LOWLYR(I,J),J) |
---|
| 215 | ! Define THSK here (for first timestep mostly) |
---|
| 216 | THSK(I,J)=TSK(I,J)/(PSFC*1.E-5)**CAPA |
---|
| 217 | ! |
---|
| 218 | !*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND) |
---|
| 219 | ! |
---|
| 220 | SEAMASK=XLAND(I,J)-1. |
---|
| 221 | ! |
---|
| 222 | !*** FILL 1-D VERTICAL ARRAYS |
---|
| 223 | !*** AND FLIP DIRECTION SINCE MYJ SCHEME |
---|
| 224 | !*** COUNTS DOWNWARD FROM THE DOMAIN'S TOP |
---|
| 225 | ! |
---|
| 226 | DO K=KTE,KTS,-1 |
---|
| 227 | KFLIP=KTE+1-K |
---|
| 228 | THK(K)=TH(I,KFLIP,J) |
---|
| 229 | TK(K)=T(I,KFLIP,J) |
---|
| 230 | RATIOMX=QV(I,KFLIP,J) |
---|
| 231 | QK(K)=RATIOMX/(1.+RATIOMX) |
---|
| 232 | PK(K)=PMID(I,KFLIP,J) |
---|
| 233 | CWMK(K)=QC(I,KFLIP,J) |
---|
| 234 | THEK(K)=(CWMK(K)*(-ELOCP/TK(K))+1.)*THK(K) |
---|
| 235 | Q2K(K)=2.*Q2(I,KFLIP,J) |
---|
| 236 | ! |
---|
| 237 | ! |
---|
| 238 | !*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES |
---|
| 239 | ! |
---|
| 240 | ZHK(K)=ZINT(I,K,J) |
---|
| 241 | ! |
---|
| 242 | ENDDO |
---|
| 243 | ZHK(KTE+1)=HT(I,J) ! Z at bottom of lowest sigma layer |
---|
| 244 | ! |
---|
| 245 | DO K=KTE,KTS,-1 |
---|
| 246 | KFLIP=KTE+1-K |
---|
| 247 | UK(K)=U(I,KFLIP,J) |
---|
| 248 | VK(K)=V(I,KFLIP,J) |
---|
| 249 | ENDDO |
---|
| 250 | ! |
---|
| 251 | !*** FIND THE HEIGHT OF THE PBL |
---|
| 252 | ! |
---|
| 253 | LPBL=LMH |
---|
| 254 | DO K=LMH-1,1,-1 |
---|
| 255 | IF(Q2K(K)<=EPSQ2L*FH) THEN |
---|
| 256 | LPBL=K |
---|
| 257 | GO TO 110 |
---|
| 258 | ENDIF |
---|
| 259 | ENDDO |
---|
| 260 | ! |
---|
| 261 | LPBL=1 |
---|
| 262 | ! |
---|
| 263 | !----------------------------------------------------------------------- |
---|
| 264 | !--------------THE HEIGHT OF THE PBL------------------------------------ |
---|
| 265 | !----------------------------------------------------------------------- |
---|
| 266 | ! |
---|
| 267 | 110 PBLH(I,J)=ZHK(LPBL)-ZHK(LMH+1) |
---|
| 268 | ! |
---|
| 269 | !---------------------------------------------------------------------- |
---|
| 270 | !*** |
---|
| 271 | !*** FIND THE SURFACE EXCHANGE COEFFICIENTS |
---|
| 272 | !*** |
---|
| 273 | !---------------------------------------------------------------------- |
---|
| 274 | PLOW=PK(LMH) |
---|
| 275 | TLOW=TK(LMH) |
---|
| 276 | THLOW=THK(LMH) |
---|
| 277 | THELOW=THEK(LMH) |
---|
| 278 | QLOW=QK(LMH) |
---|
| 279 | CWMLOW=CWMK(LMH) |
---|
| 280 | ULOW=UK(LMH) |
---|
| 281 | VLOW=VK(LMH) |
---|
| 282 | ZSL=(ZHK(LMH)-ZHK(LMH+1))*0.5 |
---|
| 283 | APESFC=(PSFC*1.E-5)**CAPA |
---|
| 284 | TZ0=THZ0(I,J)*APESFC |
---|
| 285 | ! |
---|
| 286 | CALL SFCDIF(NTSD,SEAMASK,THSK(I,J),QSFC(I,J),PSFC & |
---|
| 287 | & ,UZ0(I,J),VZ0(I,J),TZ0,THZ0(I,J),QZ0(I,J) & |
---|
| 288 | & ,USTAR(I,J),ZNT(I,J),Z0BASE(I,J),CT(I,J),RMOL(I,J) & |
---|
| 289 | & ,AKMS(I,J),AKHS(I,J),RIB(I,J),PBLH(I,J),MAVAIL(I,J) & |
---|
| 290 | & ,CHS(I,J),CHS2(I,J),CQS2(I,J) & |
---|
| 291 | & ,HFX(I,J),QFX(I,J),FLX_LH(I,J) & |
---|
| 292 | & ,FLHC(I,J),FLQC(I,J),QGH(I,J),CPM(I,J) & |
---|
| 293 | & ,ULOW,VLOW,TLOW,THLOW,THELOW,QLOW,CWMLOW & |
---|
| 294 | & ,ZSL,PLOW & |
---|
| 295 | & ,U10(I,J),V10(I,J),TSHLTR(I,J),TH10(I,J) & |
---|
| 296 | & ,QSHLTR(I,J),Q10(I,J),PSHLTR(I,J) & |
---|
| 297 | & ,IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 298 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 299 | & ,ITS,ITE,JTS,JTE,KTS,KTE,i,j,ZHK(LMH+1)) |
---|
| 300 | ! |
---|
| 301 | !*** REMOVE SUPERATURATION AT 2M AND 10M |
---|
| 302 | ! |
---|
| 303 | RAPA=APESFC |
---|
| 304 | TH02P=TSHLTR(I,J) |
---|
| 305 | TH10P=TH10(I,J) |
---|
| 306 | ! |
---|
| 307 | RAPA02=RAPA-GOCP02/TH02P |
---|
| 308 | RAPA10=RAPA-GOCP10/TH10P |
---|
| 309 | ! |
---|
| 310 | T02P=TH02P*RAPA02 |
---|
| 311 | T10P=TH10P*RAPA10 |
---|
| 312 | ! |
---|
| 313 | P02P=(RAPA02**RCAP)*1.E5 |
---|
| 314 | P10P=(RAPA10**RCAP)*1.E5 |
---|
| 315 | ! |
---|
| 316 | QS02=PQ0/P02P*EXP(A2*(T02P-A3)/(T02P-A4)) |
---|
| 317 | QS10=PQ0/P10P*EXP(A2*(T10P-A3)/(T10P-A4)) |
---|
| 318 | ! |
---|
| 319 | IF(QSHLTR(I,J)>QS02)QSHLTR(I,J)=QS02 |
---|
| 320 | IF(Q10 (I,J)>QS10)Q10 (I,J)=QS10 |
---|
| 321 | !---------------------------------------------------------------------- |
---|
| 322 | ! |
---|
| 323 | ENDDO |
---|
| 324 | ! |
---|
| 325 | !---------------------------------------------------------------------- |
---|
| 326 | ENDDO setup_integration |
---|
| 327 | !---------------------------------------------------------------------- |
---|
| 328 | |
---|
| 329 | END SUBROUTINE QNSESFC |
---|
| 330 | !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX |
---|
| 331 | !---------------------------------------------------------------------- |
---|
| 332 | SUBROUTINE SFCDIF(NTSD,SEAMASK,THS,QS,PSFC & |
---|
| 333 | & ,UZ0,VZ0,TZ0,THZ0,QZ0 & |
---|
| 334 | & ,USTAR,Z0,Z0BASE,CT,RLMO,AKMS,AKHS,RIB,PBLH,WETM & |
---|
| 335 | & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC,QGH,CPM & |
---|
| 336 | & ,ULOW,VLOW,TLOW,THLOW,THELOW,QLOW,CWMLOW & |
---|
| 337 | & ,ZSL,PLOW & |
---|
| 338 | & ,U10,V10,TH02,TH10,Q02,Q10,PSHLTR & |
---|
| 339 | & ,IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 340 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 341 | & ,ITS,ITE,JTS,JTE,KTS,KTE,i,j,ZSFC) |
---|
| 342 | ! **************************************************************** |
---|
| 343 | ! * * |
---|
| 344 | ! * SURFACE LAYER * |
---|
| 345 | ! * * |
---|
| 346 | ! **************************************************************** |
---|
| 347 | !---------------------------------------------------------------------- |
---|
| 348 | ! |
---|
| 349 | IMPLICIT NONE |
---|
| 350 | ! |
---|
| 351 | !---------------------------------------------------------------------- |
---|
| 352 | INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 353 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 354 | & ,ITS,ITE,JTS,JTE,KTS,KTE,i,j |
---|
| 355 | ! |
---|
| 356 | INTEGER,INTENT(IN) :: NTSD |
---|
| 357 | ! |
---|
| 358 | REAL,INTENT(IN) :: CWMLOW,PBLH,PLOW,QLOW,PSFC,SEAMASK,ZSFC & |
---|
| 359 | & ,THELOW,THLOW,THS,TLOW,TZ0,ULOW,VLOW,WETM,ZSL & |
---|
| 360 | & ,Z0BASE |
---|
| 361 | ! |
---|
| 362 | REAL,INTENT(OUT) :: CHS,CHS2,CPM,CQS2,CT,FLHC,FLQC,FLX_LH,HFX & |
---|
| 363 | & ,PSHLTR,Q02,Q10,QFX,QGH,RIB,RLMO,TH02,TH10,U10,V10 |
---|
| 364 | ! |
---|
| 365 | REAL,INTENT(INOUT) :: AKHS,AKMS,QZ0,THZ0,USTAR,UZ0,VZ0,Z0,QS |
---|
| 366 | !---------------------------------------------------------------------- |
---|
| 367 | !*** |
---|
| 368 | !*** LOCAL VARIABLES |
---|
| 369 | !*** |
---|
| 370 | INTEGER :: ITR,K |
---|
| 371 | ! |
---|
| 372 | REAL :: A,B,BTGH,BTGX,CXCHL,CXCHS,DTHV,DU2,ELFC,FCT & |
---|
| 373 | & ,HLFLX,HSFLX,HV,PSH02,PSH10,PSHZ,PSHZL,PSM10,PSMZ,PSMZL & |
---|
| 374 | & ,RDZ,RDZT,RLMA,RLMN,RLMP & |
---|
| 375 | & ,RLOGT,RLOGU,RWGH,RZ,RZST,RZSU,SIMH,SIMM,TEM,THM & |
---|
| 376 | & ,UMFLX,USTARK,VMFLX,WGHT,WGHTT,WGHTQ,WSTAR2 & |
---|
| 377 | & ,X,XLT,XLT4,XLU,XLU4,XT,XT4,XU,XU4,ZETALT,ZETALU & |
---|
| 378 | & ,ZETAT,ZETAU,ZQ,ZSLT,ZSLU,ZT,ZU,TOPOTERM,ZZIL |
---|
| 379 | ! |
---|
| 380 | !*** DIAGNOSTICS |
---|
| 381 | ! |
---|
| 382 | REAL :: AKHS02,AKHS10,AKMS02,AKMS10,EKMS10,QSAT10,QSAT2 & |
---|
| 383 | & ,RLNT02,RLNT10,RLNU10,SIMH02,SIMH10,SIMM10,T02,T10 & |
---|
| 384 | & ,TERM1,RLOW,U10E,V10E,WSTAR,XLT02,XLT024,XLT10 & |
---|
| 385 | & ,XLT104,XLU10,XLU104,XU10,XU104,ZT02,ZT10,ZTAT02,ZTAT10 & |
---|
| 386 | & ,ZTAU,ZTAU10,ZU10,ZUUZ |
---|
| 387 | !---------------------------------------------------------------------- |
---|
| 388 | !********************************************************************** |
---|
| 389 | !---------------------------------------------------------------------- |
---|
| 390 | RDZ=1./ZSL |
---|
| 391 | CXCHL=EXCML*RDZ |
---|
| 392 | CXCHS=EXCMS*RDZ |
---|
| 393 | ! |
---|
| 394 | BTGX=G/THLOW |
---|
| 395 | ELFC=VKARMAN*BTGX |
---|
| 396 | ! |
---|
| 397 | IF(PBLH>1000.)THEN |
---|
| 398 | BTGH=BTGX*PBLH |
---|
| 399 | ELSE |
---|
| 400 | BTGH=BTGX*1000. |
---|
| 401 | ENDIF |
---|
| 402 | ! |
---|
| 403 | !---------------------------------------------------------------------- |
---|
| 404 | ! |
---|
| 405 | !*** SEA POINTS |
---|
| 406 | ! |
---|
| 407 | !---------------------------------------------------------------------- |
---|
| 408 | ! |
---|
| 409 | IF(SEAMASK>0.5)THEN |
---|
| 410 | ! |
---|
| 411 | !---------------------------------------------------------------------- |
---|
| 412 | DO ITR=1,ITRMX |
---|
| 413 | !---------------------------------------------------------------------- |
---|
| 414 | Z0=MAX(USTFC*USTAR*USTAR,1.59E-5) |
---|
| 415 | ! |
---|
| 416 | !*** VISCOUS SUBLAYER, JANJIC MWR 1994 |
---|
| 417 | ! |
---|
| 418 | !---------------------------------------------------------------------- |
---|
| 419 | IF(USTAR<USTC)THEN |
---|
| 420 | !---------------------------------------------------------------------- |
---|
| 421 | ! |
---|
| 422 | IF(USTAR<USTR)THEN |
---|
| 423 | ! |
---|
| 424 | IF(NTSD==1)THEN |
---|
| 425 | AKMS=CXCHS |
---|
| 426 | AKHS=CXCHS |
---|
| 427 | QS=QLOW |
---|
| 428 | ENDIF |
---|
| 429 | ! |
---|
| 430 | ZU=FZU1*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR |
---|
| 431 | WGHT=AKMS*ZU*RVISC |
---|
| 432 | RWGH=WGHT/(WGHT+1.) |
---|
| 433 | UZ0=(ULOW*RWGH+UZ0)*0.5 |
---|
| 434 | VZ0=(VLOW*RWGH+VZ0)*0.5 |
---|
| 435 | ! |
---|
| 436 | ZT=FZT1*ZU |
---|
| 437 | ZQ=FZQ1*ZT |
---|
| 438 | WGHTT=AKHS*ZT*RTVISC |
---|
| 439 | WGHTQ=AKHS*ZQ*RQVISC |
---|
| 440 | ! |
---|
| 441 | IF(NTSD>1)THEN |
---|
| 442 | THZ0=((WGHTT*THLOW+THS)/(WGHTT+1.)+THZ0)*0.5 |
---|
| 443 | QZ0=((WGHTQ*QLOW+QS)/(WGHTQ+1.)+QZ0)*0.5 |
---|
| 444 | ELSE |
---|
| 445 | THZ0=(WGHTT*THLOW+THS)/(WGHTT+1.) |
---|
| 446 | QZ0=(WGHTQ*QLOW+QS)/(WGHTQ+1.) |
---|
| 447 | ENDIF |
---|
| 448 | ! |
---|
| 449 | ENDIF |
---|
| 450 | ! |
---|
| 451 | IF(USTAR>=USTR.AND.USTAR<USTC)THEN |
---|
| 452 | ZU=Z0 |
---|
| 453 | UZ0=0. |
---|
| 454 | VZ0=0. |
---|
| 455 | ! |
---|
| 456 | ZT=FZT2*SQRT(SQRT(Z0*USTAR*RVISC))/USTAR |
---|
| 457 | ZQ=FZQ2*ZT |
---|
| 458 | WGHTT=AKHS*ZT*RTVISC |
---|
| 459 | WGHTQ=AKHS*ZQ*RQVISC |
---|
| 460 | ! |
---|
| 461 | IF(NTSD>1)THEN |
---|
| 462 | THZ0=((WGHTT*THLOW+THS)/(WGHTT+1.)+THZ0)*0.5 |
---|
| 463 | QZ0=((WGHTQ*QLOW+QS)/(WGHTQ+1.)+QZ0)*0.5 |
---|
| 464 | ELSE |
---|
| 465 | THZ0=(WGHTT*THLOW+THS)/(WGHTT+1.) |
---|
| 466 | QZ0=(WGHTQ*QLOW+QS)/(WGHTQ+1.) |
---|
| 467 | ENDIF |
---|
| 468 | ! |
---|
| 469 | ENDIF |
---|
| 470 | !---------------------------------------------------------------------- |
---|
| 471 | ELSE |
---|
| 472 | !---------------------------------------------------------------------- |
---|
| 473 | ZU=Z0 |
---|
| 474 | UZ0=0. |
---|
| 475 | VZ0=0. |
---|
| 476 | ! |
---|
| 477 | ZT=Z0 |
---|
| 478 | THZ0=THS |
---|
| 479 | ! |
---|
| 480 | ZQ=Z0 |
---|
| 481 | QZ0=QS |
---|
| 482 | !---------------------------------------------------------------------- |
---|
| 483 | ENDIF |
---|
| 484 | !---------------------------------------------------------------------- |
---|
| 485 | TEM=(TLOW+TZ0)*0.5 |
---|
| 486 | THM=(THELOW+THZ0)*0.5 |
---|
| 487 | ! |
---|
| 488 | A=THM*P608 |
---|
| 489 | B=(ELOCP/TEM-1.-P608)*THM |
---|
| 490 | ! |
---|
| 491 | DTHV=((THELOW-THZ0)*((QLOW+QZ0+CWMLOW)*(0.5*P608)+1.) & |
---|
| 492 | & +(QLOW-QZ0+CWMLOW)*A+CWMLOW*B) |
---|
| 493 | ! |
---|
| 494 | DU2=MAX((ULOW-UZ0)**2+(VLOW-VZ0)**2,EPSU2) |
---|
| 495 | RIB=BTGX*DTHV*ZSL/DU2 |
---|
| 496 | !---------------------------------------------------------------------- |
---|
| 497 | ! IF(RIB>=RIC)THEN |
---|
| 498 | !---------------------------------------------------------------------- |
---|
| 499 | ! AKMS=MAX( VISC*RDZ,CXCHS) |
---|
| 500 | ! AKHS=MAX(TVISC*RDZ,CXCHS) |
---|
| 501 | !---------------------------------------------------------------------- |
---|
| 502 | ! ELSE ! turbulent branch |
---|
| 503 | !---------------------------------------------------------------------- |
---|
| 504 | ZSLU=ZSL+ZU |
---|
| 505 | ZSLT=ZSL+ZT |
---|
| 506 | ! |
---|
| 507 | RZSU=ZSLU/ZU |
---|
| 508 | RZST=ZSLT/ZT |
---|
| 509 | ! |
---|
| 510 | RLOGU=LOG(RZSU) |
---|
| 511 | RLOGT=LOG(RZST) |
---|
| 512 | ! |
---|
| 513 | !---------------------------------------------------------------------- |
---|
| 514 | !*** 1./MONIN-OBUKHOV LENGTH |
---|
| 515 | !---------------------------------------------------------------------- |
---|
| 516 | ! |
---|
| 517 | RLMO=ELFC*AKHS*DTHV/USTAR**3 |
---|
| 518 | ! |
---|
| 519 | ZETALU=ZSLU*RLMO |
---|
| 520 | ZETALT=ZSLT*RLMO |
---|
| 521 | ZETAU=ZU*RLMO |
---|
| 522 | ZETAT=ZT*RLMO |
---|
| 523 | ! |
---|
| 524 | ZETALU=MIN(MAX(ZETALU,ZTMIN1),ZTMAX1) |
---|
| 525 | ZETALT=MIN(MAX(ZETALT,ZTMIN1),ZTMAX1) |
---|
| 526 | ZETAU=MIN(MAX(ZETAU,ZTMIN1/RZSU),ZTMAX1/RZSU) |
---|
| 527 | ZETAT=MIN(MAX(ZETAT,ZTMIN1/RZST),ZTMAX1/RZST) |
---|
| 528 | ! |
---|
| 529 | !---------------------------------------------------------------------- |
---|
| 530 | !*** WATER FUNCTIONS |
---|
| 531 | !---------------------------------------------------------------------- |
---|
| 532 | ! |
---|
| 533 | RZ=(ZETAU-ZTMIN1)/DZETA1 |
---|
| 534 | K=INT(RZ) |
---|
| 535 | RDZT=RZ-REAL(K) |
---|
| 536 | K=MIN(K,KZTM2) |
---|
| 537 | K=MAX(K,0) |
---|
| 538 | PSMZ=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1) !Interpolation from the table |
---|
| 539 | ! |
---|
| 540 | RZ=(ZETALU-ZTMIN1)/DZETA1 |
---|
| 541 | K=INT(RZ) |
---|
| 542 | RDZT=RZ-REAL(K) |
---|
| 543 | K=MIN(K,KZTM2) |
---|
| 544 | K=MAX(K,0) |
---|
| 545 | PSMZL=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1) |
---|
| 546 | ! |
---|
| 547 | SIMM=PSMZL-PSMZ+RLOGU |
---|
| 548 | ! |
---|
| 549 | RZ=(ZETAT-ZTMIN1)/DZETA1 |
---|
| 550 | K=INT(RZ) |
---|
| 551 | RDZT=RZ-REAL(K) |
---|
| 552 | K=MIN(K,KZTM2) |
---|
| 553 | K=MAX(K,0) |
---|
| 554 | PSHZ=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1) |
---|
| 555 | ! |
---|
| 556 | RZ=(ZETALT-ZTMIN1)/DZETA1 |
---|
| 557 | K=INT(RZ) |
---|
| 558 | RDZT=RZ-REAL(K) |
---|
| 559 | K=MIN(K,KZTM2) |
---|
| 560 | K=MAX(K,0) |
---|
| 561 | PSHZL=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1) |
---|
| 562 | ! |
---|
| 563 | SIMH=(PSHZL-PSHZ+RLOGT)*PRT0 |
---|
| 564 | ! SIMH=(PSHZL-PSHZ+RLOGT)*FH01 |
---|
| 565 | !---------------------------------------------------------------------- |
---|
| 566 | USTARK=USTAR*VKARMAN |
---|
| 567 | AKMS=MAX(USTARK/SIMM,CXCHS) |
---|
| 568 | AKHS=MAX(USTARK/SIMH,CXCHS) |
---|
| 569 | ! |
---|
| 570 | !---------------------------------------------------------------------- |
---|
| 571 | !*** BELJAARS CORRECTION FOR USTAR (FOR CONVECTION) |
---|
| 572 | !---------------------------------------------------------------------- |
---|
| 573 | ! |
---|
| 574 | IF(DTHV<=0.)THEN !zj |
---|
| 575 | WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.) !zj |
---|
| 576 | ELSE !zj |
---|
| 577 | WSTAR2=0. !zj |
---|
| 578 | ENDIF !zj |
---|
| 579 | USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST) |
---|
| 580 | ! |
---|
| 581 | !---------------------------------------------------------------------- |
---|
| 582 | ! ENDIF ! End of turbulent branch |
---|
| 583 | !---------------------------------------------------------------------- |
---|
| 584 | ! |
---|
| 585 | ENDDO ! End of the iteration loop over sea points |
---|
| 586 | ! |
---|
| 587 | !---------------------------------------------------------------------- |
---|
| 588 | ! |
---|
| 589 | !*** LAND POINTS |
---|
| 590 | ! |
---|
| 591 | !---------------------------------------------------------------------- |
---|
| 592 | ! |
---|
| 593 | ELSE |
---|
| 594 | ! |
---|
| 595 | !---------------------------------------------------------------------- |
---|
| 596 | IF(NTSD==1)THEN |
---|
| 597 | QS=QLOW |
---|
| 598 | ENDIF |
---|
| 599 | ! |
---|
| 600 | ZU=Z0 |
---|
| 601 | UZ0=0. |
---|
| 602 | VZ0=0. |
---|
| 603 | ! |
---|
| 604 | ZT=ZU*ZTFC |
---|
| 605 | THZ0=THS |
---|
| 606 | ! |
---|
| 607 | ZQ=ZT |
---|
| 608 | QZ0=QS |
---|
| 609 | !---------------------------------------------------------------------- |
---|
| 610 | TEM=(TLOW+TZ0)*0.5 |
---|
| 611 | THM=(THELOW+THZ0)*0.5 |
---|
| 612 | ! |
---|
| 613 | A=THM*P608 |
---|
| 614 | B=(ELOCP/TEM-1.-P608)*THM |
---|
| 615 | ! |
---|
| 616 | DTHV=((THELOW-THZ0)*((QLOW+QZ0+CWMLOW)*(0.5*P608)+1.) & |
---|
| 617 | & +(QLOW-QZ0+CWMLOW)*A+CWMLOW*B) |
---|
| 618 | ! |
---|
| 619 | DU2=MAX((ULOW-UZ0)**2+(VLOW-VZ0)**2,EPSU2) |
---|
| 620 | RIB=BTGX*DTHV*ZSL/DU2 |
---|
| 621 | !---------------------------------------------------------------------- |
---|
| 622 | ! IF(RIB>=RIC)THEN |
---|
| 623 | ! AKMS=MAX( VISC*RDZ,CXCHL) |
---|
| 624 | ! AKHS=MAX(TVISC*RDZ,CXCHL) |
---|
| 625 | !---------------------------------------------------------------------- |
---|
| 626 | ! ELSE ! Turbulent branch |
---|
| 627 | !---------------------------------------------------------------------- |
---|
| 628 | ZSLU=ZSL+ZU |
---|
| 629 | ! |
---|
| 630 | RZSU=ZSLU/ZU |
---|
| 631 | ! |
---|
| 632 | RLOGU=LOG(RZSU) |
---|
| 633 | |
---|
| 634 | ZSLT=ZSL+ZU ! u,v and t are at the same level |
---|
| 635 | !---------------------------------------------------------------------- |
---|
| 636 | ! |
---|
| 637 | !mp Topo modification of ZILFC term |
---|
| 638 | ! |
---|
| 639 | TOPOTERM=TOPOFAC*ZSFC**2. |
---|
| 640 | TOPOTERM=MAX(TOPOTERM,3.0) |
---|
| 641 | ! |
---|
| 642 | IF(DTHV>0.)THEN |
---|
| 643 | ZZIL=ZILFC*TOPOTERM |
---|
| 644 | ELSE |
---|
| 645 | ZZIL=ZILFC |
---|
| 646 | ENDIF |
---|
| 647 | ! |
---|
| 648 | !---------------------------------------------------------------------- |
---|
| 649 | ! |
---|
| 650 | land_point_iteration: DO ITR=1,ITRMX |
---|
| 651 | ! |
---|
| 652 | !---------------------------------------------------------------------- |
---|
| 653 | !*** ZILITINKEVITCH FIX FOR ZT |
---|
| 654 | !---------------------------------------------------------------------- |
---|
| 655 | ! |
---|
| 656 | ! oldform ZT=MAX(EXP(ZZIL*SQRT(USTAR*ZU))*ZU,EPSZT) |
---|
| 657 | ZT=MAX(EXP(ZZIL*SQRT(USTAR*Z0BASE))*Z0BASE,EPSZT) |
---|
| 658 | ZT=MAX(ZT,ZU*ZTFC) !Quick fix by S. Sukoriansky |
---|
| 659 | ! |
---|
| 660 | RZST=ZSLT/ZT |
---|
| 661 | RLOGT=LOG(RZST) |
---|
| 662 | ! |
---|
| 663 | !---------------------------------------------------------------------- |
---|
| 664 | !*** 1./MONIN-OBUKHOV LENGTH-SCALE |
---|
| 665 | !---------------------------------------------------------------------- |
---|
| 666 | ! |
---|
| 667 | RLMO=ELFC*AKHS*DTHV/USTAR**3 |
---|
| 668 | ZETALU=ZSLU*RLMO |
---|
| 669 | ZETALT=ZSLT*RLMO |
---|
| 670 | ZETAU=ZU*RLMO |
---|
| 671 | ZETAT=ZT*RLMO |
---|
| 672 | ! |
---|
| 673 | ZETALU=MIN(MAX(ZETALU,ZTMIN2),ZTMAX2) |
---|
| 674 | ZETALT=MIN(MAX(ZETALT,ZTMIN2),ZTMAX2) |
---|
| 675 | ZETAU=MIN(MAX(ZETAU,ZTMIN2/RZSU),ZTMAX2/RZSU) |
---|
| 676 | ZETAT=MIN(MAX(ZETAT,ZTMIN2/RZST),ZTMAX2/RZST) |
---|
| 677 | ! |
---|
| 678 | !---------------------------------------------------------------------- |
---|
| 679 | !*** LAND FUNCTIONS |
---|
| 680 | !---------------------------------------------------------------------- |
---|
| 681 | ! |
---|
| 682 | RZ=(ZETAU-ZTMIN2)/DZETA2 |
---|
| 683 | K=INT(RZ) |
---|
| 684 | RDZT=RZ-REAL(K) |
---|
| 685 | K=MIN(K,KZTM2) |
---|
| 686 | K=MAX(K,0) |
---|
| 687 | PSMZ=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1) |
---|
| 688 | ! |
---|
| 689 | RZ=(ZETALU-ZTMIN2)/DZETA2 |
---|
| 690 | K=INT(RZ) |
---|
| 691 | RDZT=RZ-REAL(K) |
---|
| 692 | K=MIN(K,KZTM2) |
---|
| 693 | K=MAX(K,0) |
---|
| 694 | PSMZL=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1) |
---|
| 695 | ! |
---|
| 696 | SIMM=PSMZL-PSMZ+RLOGU |
---|
| 697 | ! |
---|
| 698 | RZ=(ZETAT-ZTMIN2)/DZETA2 |
---|
| 699 | K=INT(RZ) |
---|
| 700 | RDZT=RZ-REAL(K) |
---|
| 701 | K=MIN(K,KZTM2) |
---|
| 702 | K=MAX(K,0) |
---|
| 703 | PSHZ=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1) |
---|
| 704 | ! |
---|
| 705 | RZ=(ZETALT-ZTMIN2)/DZETA2 |
---|
| 706 | K=INT(RZ) |
---|
| 707 | RDZT=RZ-REAL(K) |
---|
| 708 | K=MIN(K,KZTM2) |
---|
| 709 | K=MAX(K,0) |
---|
| 710 | PSHZL=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1) |
---|
| 711 | ! |
---|
| 712 | SIMH=(PSHZL-PSHZ+RLOGT)*PRT0 |
---|
| 713 | ! SIMH=(PSHZL-PSHZ+RLOGT)*FH02 |
---|
| 714 | !---------------------------------------------------------------------- |
---|
| 715 | USTARK=USTAR*VKARMAN |
---|
| 716 | AKMS=MAX(USTARK/SIMM,CXCHL) |
---|
| 717 | AKHS=MAX(USTARK/SIMH,CXCHL) |
---|
| 718 | ! |
---|
| 719 | !---------------------------------------------------------------------- |
---|
| 720 | !*** BELJAARS CORRECTION FOR USTAR |
---|
| 721 | !---------------------------------------------------------------------- |
---|
| 722 | ! |
---|
| 723 | IF(DTHV<=0.)THEN !zj |
---|
| 724 | WSTAR2=WWST2*ABS(BTGH*AKHS*DTHV)**(2./3.) !zj |
---|
| 725 | ELSE !zj |
---|
| 726 | WSTAR2=0. !zj |
---|
| 727 | ENDIF !zj |
---|
| 728 | ! |
---|
| 729 | USTAR=MAX(SQRT(AKMS*SQRT(DU2+WSTAR2)),EPSUST) |
---|
| 730 | ! |
---|
| 731 | !---------------------------------------------------------------------- |
---|
| 732 | ENDDO land_point_iteration |
---|
| 733 | !---------------------------------------------------------------------- |
---|
| 734 | ! |
---|
| 735 | ! ENDIF ! End of turbulant branch over land |
---|
| 736 | ! |
---|
| 737 | !---------------------------------------------------------------------- |
---|
| 738 | ! |
---|
| 739 | ENDIF ! End of land/sea branch |
---|
| 740 | ! |
---|
| 741 | !---------------------------------------------------------------------- |
---|
| 742 | !*** COUNTERGRADIENT FIX |
---|
| 743 | !---------------------------------------------------------------------- |
---|
| 744 | ! |
---|
| 745 | ! HV=-AKHS*DTHV |
---|
| 746 | ! IF(HV>0.)THEN |
---|
| 747 | ! FCT=-10.*(BTGX)**(-1./3.) |
---|
| 748 | ! CT=FCT*(HV/(PBLH*PBLH))**(2./3.) |
---|
| 749 | ! ELSE |
---|
| 750 | CT=0. |
---|
| 751 | ! ENDIF |
---|
| 752 | ! |
---|
| 753 | !---------------------------------------------------------------------- |
---|
| 754 | !---------------------------------------------------------------------- |
---|
| 755 | !*** THE FOLLOWING DIAGNOSTIC BLOCK PRODUCES 2-m and 10-m VALUES |
---|
| 756 | !*** FOR TEMPERATURE, MOISTURE, AND WINDS. IT IS DONE HERE SINCE |
---|
| 757 | !*** THE VARIOUS QUANTITIES NEEDED FOR THE COMPUTATION ARE LOST |
---|
| 758 | !*** UPON EXIT FROM THE ROTUINE. |
---|
| 759 | !---------------------------------------------------------------------- |
---|
| 760 | !---------------------------------------------------------------------- |
---|
| 761 | ! |
---|
| 762 | WSTAR=SQRT(WSTAR2)/WWST |
---|
| 763 | ! |
---|
| 764 | UMFLX=AKMS*(ULOW -UZ0 ) |
---|
| 765 | VMFLX=AKMS*(VLOW -VZ0 ) |
---|
| 766 | HSFLX=AKHS*(THLOW-THZ0) |
---|
| 767 | HLFLX=AKHS*(QLOW -QZ0 ) |
---|
| 768 | !---------------------------------------------------------------------- |
---|
| 769 | ! IF(RIB>=RIC)THEN |
---|
| 770 | !---------------------------------------------------------------------- |
---|
| 771 | ! IF(SEAMASK>0.5)THEN |
---|
| 772 | ! AKMS10=MAX( VISC/10.,CXCHS) |
---|
| 773 | ! AKHS02=MAX(TVISC/02.,CXCHS) |
---|
| 774 | ! AKHS10=MAX(TVISC/10.,CXCHS) |
---|
| 775 | ! ELSE |
---|
| 776 | ! AKMS10=MAX( VISC/10.,CXCHL) |
---|
| 777 | ! AKHS02=MAX(TVISC/02.,CXCHL) |
---|
| 778 | ! AKHS10=MAX(TVISC/10.,CXCHL) |
---|
| 779 | ! ENDIF |
---|
| 780 | !---------------------------------------------------------------------- |
---|
| 781 | ! ELSE |
---|
| 782 | !---------------------------------------------------------------------- |
---|
| 783 | ZU10=ZU+10. |
---|
| 784 | ZT02=ZT+02. |
---|
| 785 | ZT10=ZT+10. |
---|
| 786 | ! |
---|
| 787 | RLNU10=LOG(ZU10/ZU) |
---|
| 788 | RLNT02=LOG(ZT02/ZT) |
---|
| 789 | RLNT10=LOG(ZT10/ZT) |
---|
| 790 | ! |
---|
| 791 | ZTAU10=ZU10*RLMO |
---|
| 792 | ZTAT02=ZT02*RLMO |
---|
| 793 | ZTAT10=ZT10*RLMO |
---|
| 794 | ! |
---|
| 795 | !---------------------------------------------------------------------- |
---|
| 796 | !*** SEA |
---|
| 797 | !---------------------------------------------------------------------- |
---|
| 798 | ! |
---|
| 799 | IF(SEAMASK>0.5)THEN |
---|
| 800 | ! |
---|
| 801 | !---------------------------------------------------------------------- |
---|
| 802 | ZTAU10=MIN(MAX(ZTAU10,ZTMIN1),ZTMAX1) |
---|
| 803 | ZTAT02=MIN(MAX(ZTAT02,ZTMIN1),ZTMAX1) |
---|
| 804 | ZTAT10=MIN(MAX(ZTAT10,ZTMIN1),ZTMAX1) |
---|
| 805 | !---------------------------------------------------------------------- |
---|
| 806 | RZ=(ZTAU10-ZTMIN1)/DZETA1 |
---|
| 807 | K=INT(RZ) |
---|
| 808 | RDZT=RZ-REAL(K) |
---|
| 809 | K=MIN(K,KZTM2) |
---|
| 810 | K=MAX(K,0) |
---|
| 811 | PSM10=(PSIM1(K+2)-PSIM1(K+1))*RDZT+PSIM1(K+1) |
---|
| 812 | ! |
---|
| 813 | SIMM10=PSM10-PSMZ+RLNU10 |
---|
| 814 | ! |
---|
| 815 | RZ=(ZTAT02-ZTMIN1)/DZETA1 |
---|
| 816 | K=INT(RZ) |
---|
| 817 | RDZT=RZ-REAL(K) |
---|
| 818 | K=MIN(K,KZTM2) |
---|
| 819 | K=MAX(K,0) |
---|
| 820 | PSH02=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1) |
---|
| 821 | ! |
---|
| 822 | SIMH02=(PSH02-PSHZ+RLNT02)*PRT0 |
---|
| 823 | ! SIMH02=(PSH02-PSHZ+RLNT02)*FH01 |
---|
| 824 | ! |
---|
| 825 | RZ=(ZTAT10-ZTMIN1)/DZETA1 |
---|
| 826 | K=INT(RZ) |
---|
| 827 | RDZT=RZ-REAL(K) |
---|
| 828 | K=MIN(K,KZTM2) |
---|
| 829 | K=MAX(K,0) |
---|
| 830 | PSH10=(PSIH1(K+2)-PSIH1(K+1))*RDZT+PSIH1(K+1) |
---|
| 831 | ! |
---|
| 832 | SIMH10=(PSH10-PSHZ+RLNT10)*PRT0 |
---|
| 833 | ! SIMH10=(PSH10-PSHZ+RLNT10)*FH01 |
---|
| 834 | ! |
---|
| 835 | AKMS10=MAX(USTARK/SIMM10,CXCHS) |
---|
| 836 | AKHS02=MAX(USTARK/SIMH02,CXCHS) |
---|
| 837 | AKHS10=MAX(USTARK/SIMH10,CXCHS) |
---|
| 838 | ! |
---|
| 839 | !---------------------------------------------------------------------- |
---|
| 840 | !*** LAND |
---|
| 841 | !---------------------------------------------------------------------- |
---|
| 842 | ! |
---|
| 843 | ELSE |
---|
| 844 | ! |
---|
| 845 | !---------------------------------------------------------------------- |
---|
| 846 | ZTAU10=MIN(MAX(ZTAU10,ZTMIN2),ZTMAX2) |
---|
| 847 | ZTAT02=MIN(MAX(ZTAT02,ZTMIN2),ZTMAX2) |
---|
| 848 | ZTAT10=MIN(MAX(ZTAT10,ZTMIN2),ZTMAX2) |
---|
| 849 | !---------------------------------------------------------------------- |
---|
| 850 | RZ=(ZTAU10-ZTMIN2)/DZETA2 |
---|
| 851 | K=INT(RZ) |
---|
| 852 | RDZT=RZ-REAL(K) |
---|
| 853 | K=MIN(K,KZTM2) |
---|
| 854 | K=MAX(K,0) |
---|
| 855 | PSM10=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1) |
---|
| 856 | ! |
---|
| 857 | SIMM10=PSM10-PSMZ+RLNU10 |
---|
| 858 | ! |
---|
| 859 | RZ=(ZTAT02-ZTMIN2)/DZETA2 |
---|
| 860 | K=INT(RZ) |
---|
| 861 | RDZT=RZ-REAL(K) |
---|
| 862 | K=MIN(K,KZTM2) |
---|
| 863 | K=MAX(K,0) |
---|
| 864 | PSH02=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1) |
---|
| 865 | ! |
---|
| 866 | SIMH02=(PSH02-PSHZ+RLNT02)*PRT0 |
---|
| 867 | ! SIMH02=(PSH02-PSHZ+RLNT02)*FH02 |
---|
| 868 | ! |
---|
| 869 | RZ=(ZTAT10-ZTMIN2)/DZETA2 |
---|
| 870 | K=INT(RZ) |
---|
| 871 | RDZT=RZ-REAL(K) |
---|
| 872 | K=MIN(K,KZTM2) |
---|
| 873 | K=MAX(K,0) |
---|
| 874 | PSH10=(PSIH2(K+2)-PSIH2(K+1))*RDZT+PSIH2(K+1) |
---|
| 875 | ! |
---|
| 876 | SIMH10=(PSH10-PSHZ+RLNT10)*PRT0 |
---|
| 877 | ! SIMH10=(PSH10-PSHZ+RLNT10)*FH02 |
---|
| 878 | ! |
---|
| 879 | AKMS10=MAX(USTARK/SIMM10,CXCHL) |
---|
| 880 | AKHS02=MAX(USTARK/SIMH02,CXCHL) |
---|
| 881 | AKHS10=MAX(USTARK/SIMH10,CXCHL) |
---|
| 882 | !---------------------------------------------------------------------- |
---|
| 883 | ENDIF |
---|
| 884 | !---------------------------------------------------------------------- |
---|
| 885 | ! ENDIF |
---|
| 886 | !---------------------------------------------------------------------- |
---|
| 887 | U10 =UMFLX/AKMS10+UZ0 |
---|
| 888 | V10 =VMFLX/AKMS10+VZ0 |
---|
| 889 | TH02=HSFLX/AKHS02+THZ0 |
---|
| 890 | TH10=HSFLX/AKHS10+THZ0 |
---|
| 891 | Q02 =HLFLX/AKHS02+QZ0 |
---|
| 892 | Q10 =HLFLX/AKHS10+QZ0 |
---|
| 893 | TERM1=-0.068283/TLOW |
---|
| 894 | PSHLTR=PSFC*EXP(TERM1) |
---|
| 895 | ! |
---|
| 896 | !---------------------------------------------------------------------- |
---|
| 897 | !*** COMPUTE "EQUIVALENT" Z0 TO APPROXIMATE LOCAL SHELTER READINGS. |
---|
| 898 | !---------------------------------------------------------------------- |
---|
| 899 | ! |
---|
| 900 | U10E=U10 |
---|
| 901 | V10E=V10 |
---|
| 902 | ! |
---|
| 903 | IF(SEAMASK<0.5)THEN |
---|
| 904 | |
---|
| 905 | !1st ZUUZ=MIN(0.5*ZU,0.1) |
---|
| 906 | !1st ZU=MAX(0.1*ZU,ZUUZ) |
---|
| 907 | !tst ZUUZ=amin1(ZU*0.50,0.3) |
---|
| 908 | !tst ZU=amax1(ZU*0.3,ZUUZ) |
---|
| 909 | |
---|
| 910 | ZUUZ=AMIN1(ZU*0.50,0.18) |
---|
| 911 | ZU=AMAX1(ZU*0.35,ZUUZ) |
---|
| 912 | ! |
---|
| 913 | ZU10=ZU+10. |
---|
| 914 | RZSU=ZU10/ZU |
---|
| 915 | RLNU10=LOG(RZSU) |
---|
| 916 | |
---|
| 917 | ZETAU=ZU*RLMO |
---|
| 918 | ZTAU10=ZU10*RLMO |
---|
| 919 | |
---|
| 920 | ZTAU10=MIN(MAX(ZTAU10,ZTMIN2),ZTMAX2) |
---|
| 921 | ZETAU=MIN(MAX(ZETAU,ZTMIN2/RZSU),ZTMAX2/RZSU) |
---|
| 922 | |
---|
| 923 | RZ=(ZTAU10-ZTMIN2)/DZETA2 |
---|
| 924 | K=INT(RZ) |
---|
| 925 | RDZT=RZ-REAL(K) |
---|
| 926 | K=MIN(K,KZTM2) |
---|
| 927 | K=MAX(K,0) |
---|
| 928 | PSM10=(PSIM2(K+2)-PSIM2(K+1))*RDZT+PSIM2(K+1) |
---|
| 929 | SIMM10=PSM10-PSMZ+RLNU10 |
---|
| 930 | EKMS10=MAX(USTARK/SIMM10,CXCHL) |
---|
| 931 | |
---|
| 932 | U10E=UMFLX/EKMS10+UZ0 |
---|
| 933 | V10E=VMFLX/EKMS10+VZ0 |
---|
| 934 | |
---|
| 935 | ENDIF |
---|
| 936 | ! |
---|
| 937 | U10=U10E |
---|
| 938 | V10=V10E |
---|
| 939 | ! |
---|
| 940 | !---------------------------------------------------------------------- |
---|
| 941 | !*** SET OTHER WRF DRIVER ARRAYS |
---|
| 942 | !---------------------------------------------------------------------- |
---|
| 943 | ! |
---|
| 944 | RLOW=PLOW/(R_D*TLOW) |
---|
| 945 | CHS=AKHS |
---|
| 946 | CHS2=AKHS02 |
---|
| 947 | CQS2=AKHS02 |
---|
| 948 | HFX=-RLOW*CP*HSFLX |
---|
| 949 | QFX=-RLOW*HLFLX*WETM |
---|
| 950 | FLX_LH=XLV*QFX |
---|
| 951 | FLHC=RLOW*CP*AKHS |
---|
| 952 | FLQC=RLOW*AKHS*WETM |
---|
| 953 | !!! QGH=PQ0/PSHLTR*EXP(A2S*(TSK-A3S)/(TSK-A4S)) |
---|
| 954 | QGH=((1.-SEAMASK)*PQ0+SEAMASK*PQ0SEA) & |
---|
| 955 | & /PLOW*EXP(A2S*(TLOW-A3S)/(TLOW-A4S)) |
---|
| 956 | QGH=QGH/(1.-QGH) !Convert to mixing ratio |
---|
| 957 | CPM=CP*(1.+0.8*QLOW) |
---|
| 958 | ! |
---|
| 959 | !*** DO NOT COMPUTE QS OVER LAND POINTS HERE SINCE IT IS |
---|
| 960 | !*** A PROGNOSTIC VARIABLE THERE. IT IS OKAY TO USE IT |
---|
| 961 | !*** AS A DIAGNOSTIC OVER WATER SINCE IT WILL CAUSE NO |
---|
| 962 | !*** INTERFERENCE BEFORE BEING RECOMPUTED IN MYJPBL. |
---|
| 963 | ! |
---|
| 964 | IF(SEAMASK>0.5)THEN |
---|
| 965 | QS=QLOW+QFX/(RLOW*AKHS) |
---|
| 966 | QS=QS/(1.-QS) |
---|
| 967 | ENDIF |
---|
| 968 | !---------------------------------------------------------------------- |
---|
| 969 | ! |
---|
| 970 | END SUBROUTINE SFCDIF |
---|
| 971 | ! |
---|
| 972 | !---------------------------------------------------------------------- |
---|
| 973 | SUBROUTINE QNSESFCINIT(LOWLYR,USTAR,Z0 & |
---|
| 974 | & ,SEAMASK,XICE,IVGTYP,RESTART & |
---|
| 975 | & ,ALLOWED_TO_READ & |
---|
| 976 | & ,IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 977 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 978 | & ,ITS,ITE,JTS,JTE,KTS,KTE) |
---|
| 979 | !---------------------------------------------------------------------- |
---|
| 980 | IMPLICIT NONE |
---|
| 981 | !---------------------------------------------------------------------- |
---|
| 982 | LOGICAL,INTENT(IN) :: RESTART,ALLOWED_TO_READ |
---|
| 983 | ! |
---|
| 984 | INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE & |
---|
| 985 | & ,IMS,IME,JMS,JME,KMS,KME & |
---|
| 986 | & ,ITS,ITE,JTS,JTE,KTS,KTE |
---|
| 987 | ! |
---|
| 988 | INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IVGTYP |
---|
| 989 | ! |
---|
| 990 | INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: LOWLYR |
---|
| 991 | ! |
---|
| 992 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: SEAMASK,XICE |
---|
| 993 | ! |
---|
| 994 | REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: USTAR,Z0 |
---|
| 995 | ! |
---|
| 996 | REAL,DIMENSION(0:30) :: VZ0TBL |
---|
| 997 | REAL,DIMENSION(0:30) :: VZ0TBL_24 |
---|
| 998 | ! |
---|
| 999 | INTEGER :: I,IDUM,IRECV,J,JDUM,K,ITF,JTF,KTF,MAXGBL_IVGTYP & |
---|
| 1000 | &, MAXLOC_IVGTYP,MPI_COMM_COMP |
---|
| 1001 | ! |
---|
| 1002 | ! INTEGER :: MPI_INTEGER,MPI_MAX |
---|
| 1003 | ! |
---|
| 1004 | REAL :: SM,X,ZETA1,ZETA2,ZRNG1,ZRNG2,PSIMQ,PSIHQ |
---|
| 1005 | ! |
---|
| 1006 | REAL :: PIHF=3.1415926/2.,EPS=1.E-6 |
---|
| 1007 | !---------------------------------------------------------------------- |
---|
| 1008 | VZ0TBL= & |
---|
| 1009 | & (/0., & |
---|
| 1010 | & 2.653,0.826,0.563,1.089,0.854,0.856,0.035,0.238,0.065,0.076 & |
---|
| 1011 | & ,0.011,0.035,0.011,0.000,0.000,0.000,0.000,0.000,0.000,0.000 & |
---|
| 1012 | & ,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000/) |
---|
| 1013 | |
---|
| 1014 | VZ0TBL_24= (/0., & |
---|
| 1015 | & 1.00, 0.07, 0.07, 0.07, 0.07, 0.15, & |
---|
| 1016 | & 0.08, 0.03, 0.05, 0.86, 0.80, 0.85, & |
---|
| 1017 | & 2.65, 1.09, 0.80, 0.001, 0.04, 0.05, & |
---|
| 1018 | & 0.01, 0.04, 0.06, 0.05, 0.03, 0.001, & |
---|
| 1019 | & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/) |
---|
| 1020 | |
---|
| 1021 | !---------------------------------------------------------------------- |
---|
| 1022 | ! |
---|
| 1023 | JTF=MIN0(JTE,JDE-1) |
---|
| 1024 | KTF=MIN0(KTE,KDE-1) |
---|
| 1025 | ITF=MIN0(ITE,IDE-1) |
---|
| 1026 | ! |
---|
| 1027 | ! |
---|
| 1028 | !*** FOR NOW, ASSUME SIGMA MODE FOR LOWEST MODEL LAYER |
---|
| 1029 | ! |
---|
| 1030 | DO J=JTS,JTF |
---|
| 1031 | DO I=ITS,ITF |
---|
| 1032 | LOWLYR(I,J)=1 |
---|
| 1033 | ! USTAR(I,J)=EPSUST |
---|
| 1034 | ENDDO |
---|
| 1035 | ENDDO |
---|
| 1036 | !---------------------------------------------------------------------- |
---|
| 1037 | #if (NMM_CORE == 1) |
---|
| 1038 | ! |
---|
| 1039 | IF(.NOT.RESTART)THEN |
---|
| 1040 | ! CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP) |
---|
| 1041 | ! MAXLOC_IVGTYP=MAXVAL(IVGTYP) |
---|
| 1042 | ! CALL MPI_ALLREDUCE(MAXLOC_IVGTYP,MAXGBL_IVGTYP,1,MPI_INTEGER & |
---|
| 1043 | ! &, MPI_MAX,MPI_COMM_COMP,IRECV) |
---|
| 1044 | MAXGBL_IVGTYP=MAXVAL(IVGTYP) |
---|
| 1045 | |
---|
| 1046 | #ifdef DM_PARALLEL |
---|
| 1047 | CALL WRF_DM_MAXVAL(MAXGBL_IVGTYP,IDUM,JDUM) |
---|
| 1048 | #endif |
---|
| 1049 | |
---|
| 1050 | ! |
---|
| 1051 | IF (MAXGBL_IVGTYP<13) THEN |
---|
| 1052 | DO J=JTS,JTE |
---|
| 1053 | DO I=ITS,ITE |
---|
| 1054 | SM=SEAMASK(I,J)-1. |
---|
| 1055 | IF(SM+XICE(I,J)<0.5)THEN |
---|
| 1056 | Z0(I,J)=VZ0TBL(IVGTYP(I,J)) |
---|
| 1057 | ENDIF |
---|
| 1058 | ENDDO |
---|
| 1059 | ENDDO |
---|
| 1060 | ! |
---|
| 1061 | ELSE |
---|
| 1062 | ! |
---|
| 1063 | DO J=JTS,JTE |
---|
| 1064 | DO I=ITS,ITE |
---|
| 1065 | SM=SEAMASK(I,J)-1. |
---|
| 1066 | IF(SM+XICE(I,J)<0.5)THEN |
---|
| 1067 | Z0(I,J)=VZ0TBL_24(IVGTYP(I,J)) |
---|
| 1068 | ENDIF |
---|
| 1069 | ENDDO |
---|
| 1070 | ENDDO |
---|
| 1071 | ! |
---|
| 1072 | ENDIF ! Vegtype check |
---|
| 1073 | ! |
---|
| 1074 | ENDIF ! Restart check |
---|
| 1075 | #endif |
---|
| 1076 | !---------------------------------------------------------------------- |
---|
| 1077 | IF(.NOT.RESTART)THEN |
---|
| 1078 | DO J=JTS,JTE |
---|
| 1079 | DO I=ITS,ITF |
---|
| 1080 | USTAR(I,J)=0.1 |
---|
| 1081 | ENDDO |
---|
| 1082 | ENDDO |
---|
| 1083 | ENDIF |
---|
| 1084 | !---------------------------------------------------------------------- |
---|
| 1085 | ! |
---|
| 1086 | !*** COMPUTE SURFACE LAYER INTEGRAL FUNCTIONS |
---|
| 1087 | ! |
---|
| 1088 | !---------------------------------------------------------------------- |
---|
| 1089 | FH01=1. |
---|
| 1090 | FH02=1. |
---|
| 1091 | ! |
---|
| 1092 | ! ZTMIN1=-10.0 |
---|
| 1093 | ! ZTMAX1=2.0 |
---|
| 1094 | ! ZTMIN2=-10.0 |
---|
| 1095 | ! ZTMAX2=2.0 |
---|
| 1096 | ZTMIN1=-5.0 |
---|
| 1097 | ! ZTMAX1=3.0 |
---|
| 1098 | ZTMAX1=2.3 ! changed |
---|
| 1099 | ZTMIN2=-5.0 |
---|
| 1100 | ! ZTMAX2=3.0 |
---|
| 1101 | ZTMAX2=2.3 ! changed |
---|
| 1102 | ! |
---|
| 1103 | ZRNG1=ZTMAX1-ZTMIN1 |
---|
| 1104 | ZRNG2=ZTMAX2-ZTMIN2 |
---|
| 1105 | ! |
---|
| 1106 | DZETA1=ZRNG1/(KZTM-1) |
---|
| 1107 | DZETA2=ZRNG2/(KZTM-1) |
---|
| 1108 | ! |
---|
| 1109 | !---------------------------------------------------------------------- |
---|
| 1110 | !*** FUNCTION DEFINITION LOOP |
---|
| 1111 | !---------------------------------------------------------------------- |
---|
| 1112 | ! |
---|
| 1113 | ZETA1=ZTMIN1 |
---|
| 1114 | ZETA2=ZTMIN2 |
---|
| 1115 | ! |
---|
| 1116 | DO K=1,KZTM |
---|
| 1117 | ! |
---|
| 1118 | !---------------------------------------------------------------------- |
---|
| 1119 | !*** UNSTABLE RANGE |
---|
| 1120 | !---------------------------------------------------------------------- |
---|
| 1121 | ! |
---|
| 1122 | IF(ZETA1<0.)THEN |
---|
| 1123 | ! |
---|
| 1124 | !---------------------------------------------------------------------- |
---|
| 1125 | !*** PAULSON 1970 FUNCTIONS |
---|
| 1126 | !---------------------------------------------------------------------- |
---|
| 1127 | X=SQRT(SQRT(1.-16.*ZETA1)) |
---|
| 1128 | |
---|
| 1129 | PSIM1(K)=-2.*LOG((X+1.)/2.)-LOG((X*X+1.)/2.)+2.*ATAN(X)-PIHF |
---|
| 1130 | PSIH1(K)=-2.*LOG((X*X+1.)/2.) |
---|
| 1131 | ! |
---|
| 1132 | !---------------------------------------------------------------------- |
---|
| 1133 | !*** DERIVED BY S. SUKORIANSKY USING QNSE THEORY |
---|
| 1134 | !---------------------------------------------------------------------- |
---|
| 1135 | ! |
---|
| 1136 | PSIMQ=2.25*ZETA1+ZETA1**2+45.9*ZETA1**3 |
---|
| 1137 | PSIHQ=2.0*ZETA1+10.3*ZETA1**2+130*ZETA1**3 |
---|
| 1138 | ! PSIHQ=1.45*ZETA1+7.42*ZETA1**2+93.5*ZETA1**3 ! /PRT0 |
---|
| 1139 | PSIM1(K)=MAX(PSIMQ,PSIM1(K)) |
---|
| 1140 | PSIH1(K)=MAX(PSIHQ,PSIH1(K)) |
---|
| 1141 | ! |
---|
| 1142 | !---------------------------------------------------------------------- |
---|
| 1143 | !*** STABLE RANGE |
---|
| 1144 | !---------------------------------------------------------------------- |
---|
| 1145 | ! |
---|
| 1146 | ELSE |
---|
| 1147 | ! |
---|
| 1148 | !---------------------------------------------------------------------- |
---|
| 1149 | !*** PAULSON 1970 FUNCTIONS |
---|
| 1150 | !---------------------------------------------------------------------- |
---|
| 1151 | ! |
---|
| 1152 | ! PSIM1(K)=5.*ZETA1 |
---|
| 1153 | ! PSIH1(K)=5.*ZETA1 |
---|
| 1154 | !---------------------------------------------------------------------- |
---|
| 1155 | !*** HOLTSLAG AND DE BRUIN 1988 |
---|
| 1156 | !---------------------------------------------------------------------- |
---|
| 1157 | ! |
---|
| 1158 | ! PSIM1(K)=0.7*ZETA1+0.75*ZETA1*(6.-0.35*ZETA1)*EXP(-0.35*ZETA1) |
---|
| 1159 | ! PSIH1(K)=0.7*ZETA1+0.75*ZETA1*(6.-0.35*ZETA1)*EXP(-0.35*ZETA1) |
---|
| 1160 | ! |
---|
| 1161 | !---------------------------------------------------------------------- |
---|
| 1162 | !*** DERIVED BY S. SUKORIANSKY USING QNSE THEORY |
---|
| 1163 | !---------------------------------------------------------------------- |
---|
| 1164 | ! |
---|
| 1165 | PSIM1(K)=2.25*ZETA1-0.2*ZETA1*ZETA1 |
---|
| 1166 | PSIH1(K)=2.0*ZETA1+0.14*((ZETA1-0.5)**5+0.03125) |
---|
| 1167 | ! PSIH1(K)=1.42*ZETA1+0.1*((ZETA1-0.5)**5+0.03125) !/PRT0 |
---|
| 1168 | !---------------------------------------------------------------------- |
---|
| 1169 | ! |
---|
| 1170 | ENDIF |
---|
| 1171 | ! |
---|
| 1172 | !---------------------------------------------------------------------- |
---|
| 1173 | !*** UNSTABLE RANGE |
---|
| 1174 | !---------------------------------------------------------------------- |
---|
| 1175 | ! |
---|
| 1176 | IF(ZETA2<0.)THEN |
---|
| 1177 | ! |
---|
| 1178 | !---------------------------------------------------------------------- |
---|
| 1179 | !*** PAULSON 1970 FUNCTIONS |
---|
| 1180 | !---------------------------------------------------------------------- |
---|
| 1181 | ! |
---|
| 1182 | X=SQRT(SQRT(1.-16.*ZETA2)) |
---|
| 1183 | |
---|
| 1184 | PSIM2(K)=-2.*LOG((X+1.)/2.)-LOG((X*X+1.)/2.)+2.*ATAN(X)-PIHF |
---|
| 1185 | PSIH2(K)=-2.*LOG((X*X+1.)/2.) |
---|
| 1186 | ! |
---|
| 1187 | !---------------------------------------------------------------------- |
---|
| 1188 | !*** DERIVED BY S. SUKORIANSKY USING QNSE THEORY |
---|
| 1189 | !---------------------------------------------------------------------- |
---|
| 1190 | ! |
---|
| 1191 | PSIMQ=2.25*ZETA2+ZETA2**2+45.9*ZETA2**3 |
---|
| 1192 | PSIHQ=2.0*ZETA2+10.3*ZETA2**2+130*ZETA2**3 |
---|
| 1193 | ! PSIHQ=1.45*ZETA2+7.42*ZETA2**2+93.5*ZETA2**3 ! /PRT0 |
---|
| 1194 | PSIM2(K)=MAX(PSIMQ,PSIM2(K)) |
---|
| 1195 | PSIH2(K)=MAX(PSIHQ,PSIH2(K)) |
---|
| 1196 | ! |
---|
| 1197 | !---------------------------------------------------------------------- |
---|
| 1198 | !*** STABLE RANGE |
---|
| 1199 | !---------------------------------------------------------------------- |
---|
| 1200 | ! |
---|
| 1201 | ELSE |
---|
| 1202 | ! |
---|
| 1203 | !---------------------------------------------------------------------- |
---|
| 1204 | !*** PAULSON 1970 FUNCTIONS |
---|
| 1205 | !---------------------------------------------------------------------- |
---|
| 1206 | ! |
---|
| 1207 | ! PSIM2(K)=5.*ZETA2 |
---|
| 1208 | ! PSIH2(K)=5.*ZETA2 |
---|
| 1209 | ! |
---|
| 1210 | !---------------------------------------------------------------------- |
---|
| 1211 | !*** HOLTSLAG AND DE BRUIN 1988 |
---|
| 1212 | !---------------------------------------------------------------------- |
---|
| 1213 | ! |
---|
| 1214 | ! PSIM2(K)=0.7*ZETA2+0.75*ZETA2*(6.-0.35*ZETA2)*EXP(-0.35*ZETA2) |
---|
| 1215 | ! PSIH2(K)=0.7*ZETA2+0.75*ZETA2*(6.-0.35*ZETA2)*EXP(-0.35*ZETA2) |
---|
| 1216 | ! |
---|
| 1217 | !---------------------------------------------------------------------- |
---|
| 1218 | !*** DERIVED BY S. SUKORIANSKY USING QNSE THEORY |
---|
| 1219 | !---------------------------------------------------------------------- |
---|
| 1220 | ! |
---|
| 1221 | PSIM2(K)=2.25*ZETA2-0.2*ZETA2*ZETA2 |
---|
| 1222 | PSIH2(K)=2.0*ZETA2+0.14*((ZETA2-0.5)**5+0.03125) |
---|
| 1223 | ! PSIH2(K)=1.42*ZETA2+0.1*((ZETA2-0.5)**5+0.03125) !/PRT0 |
---|
| 1224 | !---------------------------------------------------------------------- |
---|
| 1225 | ! |
---|
| 1226 | ENDIF |
---|
| 1227 | ! |
---|
| 1228 | !---------------------------------------------------------------------- |
---|
| 1229 | IF(K==KZTM)THEN |
---|
| 1230 | ZTMAX1=ZETA1 |
---|
| 1231 | ZTMAX2=ZETA2 |
---|
| 1232 | ENDIF |
---|
| 1233 | ! |
---|
| 1234 | ZETA1=ZETA1+DZETA1 |
---|
| 1235 | ZETA2=ZETA2+DZETA2 |
---|
| 1236 | !---------------------------------------------------------------------- |
---|
| 1237 | ENDDO |
---|
| 1238 | !---------------------------------------------------------------------- |
---|
| 1239 | ZTMAX1=ZTMAX1-EPS |
---|
| 1240 | ZTMAX2=ZTMAX2-EPS |
---|
| 1241 | !---------------------------------------------------------------------- |
---|
| 1242 | ! |
---|
| 1243 | END SUBROUTINE QNSESFCINIT |
---|
| 1244 | ! |
---|
| 1245 | !---------------------------------------------------------------------- |
---|
| 1246 | ! |
---|
| 1247 | END MODULE MODULE_SF_QNSESFC |
---|
| 1248 | ! |
---|
| 1249 | !---------------------------------------------------------------------- |
---|