Changeset 1377 for trunk/LMDZ.MARS/libf
- Timestamp:
- Feb 13, 2015, 8:58:27 AM (10 years ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 2 edited
-
pbl_parameters.F (modified) (7 diffs)
-
physiq.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F
r1226 r1377 145 145 ite=0. 146 146 residual=abs(pz0tcomp(ig)-pz0t) 147 ! zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 148 ! & ,(0.3*wstar_in(ig))**2) 147 149 148 zu2(ig)=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 150 149 & + (log(1.+0.7*wstar_in(ig) + 2.3*wstar_in(ig)**2))**2 … … 158 157 & *sqrt(zzlay(ig,1)*pz0(ig)) 159 158 & *(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t))) 160 & *(ph(ig,1)-pts(ig))/ zu2(ig)159 & *(ph(ig,1)-pts(ig))/(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)) 161 160 ELSE 162 161 print*,'warning, infinite Richardson at surface' … … 220 219 221 220 222 ! u* theta* computation 223 ! and Monin Obukhov length: 221 ! u* theta* computation: 224 222 225 223 DO ig=1,ngrid … … 227 225 ustar(ig)=0. 228 226 tstar(ig)=0. 229 L_mo(ig)=0. 227 230 228 ELSE 231 229 ustar(ig)=sqrt(pcdv(ig)) … … 233 231 tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1)) 234 232 & /sqrt(pcdv(ig)) 235 L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL 236 ENDIF 237 ENDDO 238 239 ! Monin Obukhov length: 240 241 ! DO ig=1,ngrid 242 ! IF (rib(ig) .ge. ric) THEN 243 ! L_mo(ig)=0. 244 ! ELSE 245 ! L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL 246 ! ENDIF 247 ! ENDDO 248 249 ! Interpolation: 233 234 ENDIF 235 ENDDO 250 236 251 237 DO ig=1,ngrid … … 253 239 u_out(ig,n)=0. 254 240 Teta_out(ig,n)=pts(ig) 255 ELSEIF (L_mo(ig) .gt. 0.) THEN 256 ! u_out(ig,n)=(ustar(ig)/karman)*log(zout/pz0(ig)) + 257 ! & 5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig)) 258 ! Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) 259 ! & *log(zout/pz0tcomp(ig)) + 260 ! & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) 261 ! & *(zout-pz0tcomp(ig)) 262 IF ((zout/L_mo(ig).lt.ric).and.(pz0(ig).lt.ric)) THEN 263 ! & ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) ) THEN 264 265 u_out(ig,n)=(ustar(ig)/karman)*( 266 & log((ric-pz0(ig)/L_mo(ig))/(ric-zout/L_mo(ig))) 267 & +log(zout/pz0(ig)) 268 & ) 269 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) THEN 270 271 u_out(ig,n)=(ustar(ig)/karman)*( 272 & log((zout-ric*L_mo(ig))/(pz0(ig)-ric*L_mo(ig))) 273 & +log(pz0(ig)/zout) 274 & ) 275 276 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).lt.ric)) THEN 277 !integral of the stability function does not converge 278 279 280 u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) 281 282 283 ENDIF 284 IF((zout/L_mo(ig).lt.ric).and.(pz0tcomp(ig).lt.ric)) THEN 285 ! & ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) ) THEN 286 287 Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( 288 & log((ric-pz0tcomp(ig)/L_mo(ig))/(ric-zout/L_mo(ig))) 289 & +log(zout/pz0tcomp(ig)) 290 & ) 291 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) THEN 292 293 Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( 294 & log((zout-ric*L_mo(ig))/(pz0tcomp(ig)-ric*L_mo(ig))) 295 & +log(pz0tcomp(ig)/zout) 296 & ) 297 298 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).lt.ric)) THEN 299 !integral of the stability function does not converge 300 301 Teta_out(ig,n)=ph(ig,1) 302 303 ENDIF 304 305 ELSEIF (L_mo(ig) .lt. 0.) THEN 306 307 IF(L_mo(ig) .gt. -1000.) THEN 308 309 u_out(ig,n)=(ustar(ig)/karman)*( 310 & 2.*atan((1.-16.*zout/L_mo(ig))**0.25) 311 & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25) 312 & -2.*log(1.+(1.-16.*zout/L_mo(ig))**0.25) 313 & +2.*log(1.+(1.-16.*pz0(ig)/L_mo(ig))**0.25) 314 & - log(1.+sqrt(1.-16.*zout/L_mo(ig))) 315 & + log(1.+sqrt(1.-16.*pz0(ig)/L_mo(ig))) 316 & + log(zout/pz0(ig)) 317 & ) 318 319 Teta_out(ig,n)=MAX(pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( 320 & 2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig))) 321 & -2.*log(1.+sqrt(1.-16.*zout/L_mo(ig))) 322 & + log(zout/pz0tcomp(ig)) 323 & ),ph(ig,1)) 324 325 ELSE 326 327 ! We have to treat the case where L is very large and negative, 328 ! corresponding to a very nearly stable atmosphere (but not quite !) 329 ! i.e. teta* <0 and almost zero. We choose the cutoff 330 ! corresponding to L_mo < -1000 and use a 3rd order taylor expansion. The difference 331 ! between the two expression for z/L = -1/1000 is -1.54324*10^-9 332 ! (we do that to avoid using r*4 precision, otherwise, we get -inf values) 333 334 u_out(ig,n)=(ustar(ig)/karman)*( 335 & (4./L_mo(ig))*(zout-pz0(ig)) 336 & + (20./(L_mo(ig))**2)*(zout**2-pz0(ig)**2) 337 & + (160./(L_mo(ig))**3)*(zout**3-pz0(ig)**3) 338 & + log(zout/pz0(ig)) 339 & ) 340 341 Teta_out(ig,n)=MAX(pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( 342 & (8./L_mo(ig))*(zout-pz0tcomp(ig)) 343 & + (48./(L_mo(ig))**2)*(zout**2-pz0tcomp(ig)**2) 344 & + (1280./(3.*(L_mo(ig))**3))*(zout**3-pz0tcomp(ig)**3) 345 & + log(zout/pz0tcomp(ig)) 346 & ),ph(ig,1)) 347 348 ENDIF 349 ELSE 350 u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) 351 Teta_out(ig,n)=ph(ig,1) 241 ELSE 242 u_out(ig,n)= ustar(ig)*log(zout/pz0(ig))/ 243 &(karman*sqrt(fm(ig))) 244 245 Teta_out(ig,n)=pts(ig)+(tstar(ig)*sqrt(fm(ig))*log(zout/ 246 & (pz0tcomp(ig)))/ 247 &(karman*fh(ig))) 248 352 249 ENDIF 353 IF(zout .lt. pz0(ig)) THEN 250 251 IF (zout .lt. pz0(ig)) THEN 354 252 u_out(ig,n)=0. 355 ENDIF 356 357 ! Final check 358 IF (L_mo(ig) .gt. 0) THEN 359 IF (Teta_out(ig,n) .gt. ph(ig,1)) THEN 360 Teta_out(ig,n)=ph(ig,1) 361 ELSEIF (Teta_out(ig,n) .lt. pts(ig)) THEN 362 Teta_out(ig,n)=pts(ig) 363 ENDIF 364 ELSEIF (L_mo(ig) .lt. 0) THEN 365 IF (Teta_out(ig,n) .lt. ph(ig,1)) THEN 366 Teta_out(ig,n)=ph(ig,1) 367 ELSEIF (Teta_out(ig,n) .gt. pts(ig)) THEN 368 Teta_out(ig,n)=pts(ig) 369 ENDIF 370 ENDIF 371 372 IF (u_out(ig,n) .gt. sqrt(pu(ig,1)**2 + pv(ig,1)**2)) THEN 373 u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) 374 ENDIF 253 ENDIF 375 254 376 255 ENDDO … … 488 367 & 'momentum stab function','m/s',2,fh) 489 368 call WRITEDIAGFI(ngrid,'B_pbl', 490 & ' flottabilité','m/',2,(ph(:,1)-pts(:))/pts(:))369 & 'buoyancy','m/',2,(ph(:,1)-pts(:))/pts(:)) 491 370 call WRITEDIAGFI(ngrid,'zot_pbl', 492 & ' flottabilité','ms',2,pz0tcomp)493 call WRITEDIAGFI(ngrid,'zz1',' flottabilité','m',2,zzlay(:,1))371 & 'buoyancy','ms',2,pz0tcomp) 372 call WRITEDIAGFI(ngrid,'zz1','buoyancy','m',2,zzlay(:,1)) 494 373 #endif 495 374 -
trunk/LMDZ.MARS/libf/phymars/physiq.F
r1375 r1377 518 518 zdtsurf(:)=0 519 519 dqsurf(:,:)=0 520 #ifdef DUSTSTORM 520 521 pq_tmp(:,:,:)=0 522 #endif 521 523 igout=ngrid/2+1 522 524 … … 1567 1569 call WRITEDIAGFI(ngrid,'teta_star', 1568 1570 & 'friction potential temperature','K',dimout,tstar) 1569 call WRITEDIAGFI(ngrid,'L',1570 & 'Monin Obukhov length','m',dimout,L_mo)1571 ! call WRITEDIAGFI(ngrid,'L', 1572 ! & 'Monin Obukhov length','m',dimout,L_mo) 1571 1573 call WRITEDIAGFI(ngrid,'vvv', 1572 1574 & 'Vertical velocity variance at zout','m',dimout,vvv)
Note: See TracChangeset
for help on using the changeset viewer.
