Changeset 636
- Timestamp:
- Apr 27, 2012, 11:22:46 AM (13 years ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F
r605 r636 144 144 lambda=(sqrt(bh/bm))/alphah 145 145 ric=betah/(betam**2) 146 147 146 DO ig=1,ngrid 148 147 ite=0. 149 148 residual=abs(pz0tcomp(ig)-pz0t) 150 zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 151 & ,(0.3*wstar_in(ig))**2) 149 ! zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 150 ! & ,(0.3*wstar_in(ig))**2) 151 zu2(ig)=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)+ 152 & + (log(1.+0.7*wstar_in(ig) + 2.3*wstar_in(ig)**2))**2 152 153 153 154 DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.)) … … 157 158 ! Richardson number formulation proposed by D.E. England et al. (1995) 158 159 rib(ig) = (pg/pts(ig)) 159 & *sqrt(zzl ev(ig,2)*pz0(ig))160 & *(((log(zzl ev(ig,2)/pz0(ig)))**2)/(log(zzlev(ig,2)/pz0t)))160 & *sqrt(zzlay(ig,1)*pz0(ig)) 161 & *(((log(zzlay(ig,1)/pz0(ig)))**2)/(log(zzlay(ig,1)/pz0t))) 161 162 & *(ph(ig,1)-pts(ig))/zu2(ig) 162 163 ELSE … … 166 167 ENDIF 167 168 168 z1z0=zzl ev(ig,2)/pz0(ig)169 z1z0t=zzl ev(ig,2)/pz0t169 z1z0=zzlay(ig,1)/pz0(ig) 170 z1z0t=zzlay(ig,1)/pz0t 170 171 171 172 cdn(ig)=karman/log(z1z0) … … 221 222 222 223 223 ! Large-scale wind at first layer (without gustiness) and224 224 ! u* theta* computation 225 ! and Monin Obukhov length: 225 226 226 227 DO ig=1,ngrid … … 228 229 ustar(ig)=0. 229 230 tstar(ig)=0. 231 L_mo(ig)=0. 230 232 ELSE 231 233 ustar(ig)=sqrt(pcdv(ig)) … … 233 235 tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1)) 234 236 & /sqrt(pcdv(ig)) 237 L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL 235 238 ENDIF 236 239 ENDDO … … 238 241 ! Monin Obukhov length: 239 242 240 DO ig=1,ngrid241 IF (rib(ig) .gt. ric) THEN242 L_mo(ig)=0.243 ELSE244 L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL245 ENDIF246 ENDDO243 ! DO ig=1,ngrid 244 ! IF (rib(ig) .ge. ric) THEN 245 ! L_mo(ig)=0. 246 ! ELSE 247 ! L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig)) !as defined here, L is positive for SBL, negative for UBL 248 ! ENDIF 249 ! ENDDO 247 250 248 251 ! Interpolation: … … 253 256 Teta_out(ig,n)=pts(ig) 254 257 ELSEIF (L_mo(ig) .gt. 0.) THEN 255 u_out(ig,n)=(ustar(ig)/karman)*log(zout/pz0(ig)) + 256 & 5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig)) 257 Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) 258 & *log(zout/pz0tcomp(ig)) + 259 & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) 260 & *(zout-pz0tcomp(ig)) 258 ! u_out(ig,n)=(ustar(ig)/karman)*log(zout/pz0(ig)) + 259 ! & 5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig)) 260 ! Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) 261 ! & *log(zout/pz0tcomp(ig)) + 262 ! & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) 263 ! & *(zout-pz0tcomp(ig)) 264 IF ((zout/L_mo(ig).lt.ric).and.(pz0(ig).lt.ric)) THEN 265 ! & ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) ) THEN 266 267 u_out(ig,n)=(ustar(ig)/karman)*( 268 & log((ric-pz0(ig)/L_mo(ig))/(ric-zout/L_mo(ig))) 269 & +log(zout/pz0(ig)) 270 & ) 271 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).gt.ric)) THEN 272 273 u_out(ig,n)=(ustar(ig)/karman)*( 274 & log((zout-ric*L_mo(ig))/(pz0(ig)-ric*L_mo(ig))) 275 & +log(pz0(ig)/zout) 276 & ) 277 278 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0(ig).lt.ric)) THEN 279 !integral of the stability function does not converge 280 281 282 u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) 283 284 285 ENDIF 286 IF((zout/L_mo(ig).lt.ric).and.(pz0tcomp(ig).lt.ric)) THEN 287 ! & ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) ) THEN 288 289 Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( 290 & log((ric-pz0tcomp(ig)/L_mo(ig))/(ric-zout/L_mo(ig))) 291 & +log(zout/pz0tcomp(ig)) 292 & ) 293 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).gt.ric)) THEN 294 295 Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( 296 & log((zout-ric*L_mo(ig))/(pz0tcomp(ig)-ric*L_mo(ig))) 297 & +log(pz0tcomp(ig)/zout) 298 & ) 299 300 ELSEIF ((zout/L_mo(ig).gt.ric).and.(pz0tcomp(ig).lt.ric)) THEN 301 !integral of the stability function does not converge 302 303 Teta_out(ig,n)=ph(ig,1) 304 305 ENDIF 306 261 307 ELSEIF (L_mo(ig) .lt. 0.) THEN 262 308 … … 304 350 ENDIF 305 351 ELSE 306 u_out(ig,n)= 0.307 Teta_out(ig,n)=p ts(ig)352 u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) 353 Teta_out(ig,n)=ph(ig,1) 308 354 ENDIF 309 355 IF(zout .lt. pz0(ig)) THEN 310 356 u_out(ig,n)=0. 311 357 ENDIF 358 359 ! Final check 360 IF (L_mo(ig) .gt. 0) THEN 361 IF (Teta_out(ig,n) .gt. ph(ig,1)) THEN 362 Teta_out(ig,n)=ph(ig,1) 363 ELSEIF (Teta_out(ig,n) .lt. pts(ig)) THEN 364 Teta_out(ig,n)=pts(ig) 365 ENDIF 366 ELSEIF (L_mo(ig) .lt. 0) THEN 367 IF (Teta_out(ig,n) .lt. ph(ig,1)) THEN 368 Teta_out(ig,n)=ph(ig,1) 369 ELSEIF (Teta_out(ig,n) .gt. pts(ig)) THEN 370 Teta_out(ig,n)=pts(ig) 371 ENDIF 372 ENDIF 373 374 IF (u_out(ig,n) .gt. sqrt(pu(ig,1)**2 + pv(ig,1)**2)) THEN 375 u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2) 376 ENDIF 377 312 378 ENDDO 313 379 … … 420 486 ENDIF ! of if calltherm 421 487 488 call WRITEDIAGFI(ngrid,'rib_pbl', 489 & 'Richardson in pbl parameter','m/s',2,rib) 490 call WRITEDIAGFI(ngrid,'cdn_pbl', 491 & 'neutral momentum coef','m/s',2,cdn) 492 call WRITEDIAGFI(ngrid,'fm_pbl', 493 & 'momentum stab function','m/s',2,fm) 494 call WRITEDIAGFI(ngrid,'uv', 495 & 'wind norm first layer','m/s',2,sqrt(zu2)) 496 call WRITEDIAGFI(ngrid,'uvtrue', 497 & 'wind norm first layer','m/s',2,sqrt(pu(:,1)**2+pv(:,1)**2)) 498 call WRITEDIAGFI(ngrid,'chn_pbl', 499 & 'neutral momentum coef','m/s',2,chn) 500 call WRITEDIAGFI(ngrid,'fh_pbl', 501 & 'momentum stab function','m/s',2,fh) 502 call WRITEDIAGFI(ngrid,'B_pbl', 503 & 'flottabilité','m/',2,(ph(:,1)-pts(:))/pts(:)) 504 call WRITEDIAGFI(ngrid,'zot_pbl', 505 & 'flottabilité','ms',2,pz0tcomp) 506 call WRITEDIAGFI(ngrid,'zz1','flottabilité','m',2,zzlay(:,1)) 507 422 508 RETURN 423 509 END -
trunk/LMDZ.MARS/libf/phymars/physiq.F
r635 r636 347 347 REAL dtke_th(ngridmx,nlayermx+1) 348 348 REAL zcdv(ngridmx), zcdh(ngridmx) 349 REAL, ALLOCATABLE, DIMENSION(:,:) :: T eta_out349 REAL, ALLOCATABLE, DIMENSION(:,:) :: T_out 350 350 REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out 351 REAL u_out1(ngridmx), T_out1(ngridmx) 351 352 REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters] 352 353 REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp … … 354 355 REAL zu2(ngridmx) 355 356 c======================================================================= 356 357 print*,'physiq0', nq,nqmx358 359 357 360 358 c 1. Initialisation: … … 1364 1362 c ---------------------------------------------------------- 1365 1363 1366 IF (1 .eq. 0.) THEN 1367 IF (callrichsl) THEN 1368 n_out=5 1369 1370 ALLOCATE(z_out(n_out)) 1371 ALLOCATE(Teta_out(ngrid,n_out)) 1372 ALLOCATE(u_out(ngrid,n_out)) 1373 1374 z_out(:)=[0.001,0.05,0.1,0.5,1.] 1375 u_out(:,:)=0. 1376 Teta_out(:,:)=0. 1377 1378 call pbl_parameters(ngrid,nlayer,ps,zplay,z0, 1364 n_out=5 ! number of elements in the z_out array. 1365 ! for z_out=[3.,2.,1.,0.5,0.1], n_out must be set 1366 ! to 5 1367 IF (n_out .ne. 0) THEN 1368 1369 IF(.NOT. ALLOCATED(z_out)) ALLOCATE(z_out(n_out)) 1370 IF(.NOT. ALLOCATED(T_out)) ALLOCATE(T_out(ngrid,n_out)) 1371 IF(.NOT. ALLOCATED(u_out)) ALLOCATE(u_out(ngrid,n_out)) 1372 1373 z_out(:)=[3.,2.,1.,0.5,0.1] 1374 u_out(:,:)=0. 1375 T_out(:,:)=0. 1376 1377 call pbl_parameters(ngrid,nlayer,ps,zplay,z0, 1379 1378 & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out, 1380 & T eta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv)1379 & T_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) 1381 1380 1382 1381 #ifndef MESOSCALE … … 1387 1386 ENDIF 1388 1387 DO n=1,n_out 1389 write(zstring, '(F 9.6)') z_out(n)1390 call WRITEDIAGFI(ngrid,'T eta_out_'//trim(zstring),1391 & 'potential temperature at z_out','K',dimout,T eta_out(:,n))1388 write(zstring, '(F8.6)') z_out(n) 1389 call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring), 1390 & 'potential temperature at z_out','K',dimout,T_out(:,n)) 1392 1391 call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring), 1393 1392 & 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n)) 1394 ENDDO 1393 ENDDO 1395 1394 call WRITEDIAGFI(ngrid,'u_star', 1396 1395 & 'friction velocity','m/s',dimout,ustar) … … 1403 1402 call WRITEDIAGFI(ngrid,'vhf', 1404 1403 & 'Vertical heat flux at zout','m',dimout,vhf) 1404 #else 1405 T_out1(:)=T_out(:,1) 1406 u_out1(:)=u_out(:,1) 1405 1407 #endif 1406 1408 1407 1409 ENDIF 1408 ENDIF ! of pbl interpolation outputs1409 1410 1410 1411 c ---------------------------------------------------------- -
trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
r604 r636 148 148 ! zu2=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1),0.25*wstar(ig)**2) 149 149 ! zu2=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) 150 zu2(ig)=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1), & 151 & (0.3*wstar(ig))**2) 150 ! zu2(ig)=MAX(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1), & 151 ! & (0.3*wstar(ig))**2) 152 zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + 153 & + (log(1.+0.7*wstar(ig) + 2.3*wstar(ig)**2))**2 154 152 155 ! zu2(ig)=pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (0.5*wstar(ig))**2 153 156
Note: See TracChangeset
for help on using the changeset viewer.
