- Timestamp:
- Feb 14, 2012, 12:41:24 PM (13 years ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F
r496 r529 1 SUBROUTINE pbl_parameters(ngrid,nlay,p z0,2 & pg, pz,pu,pv,wmax,hfmax,zmax,pts,ph,z_out,3 & Teta_out,u_out,ustar,tstar, wstar,L_mo,vhf,vvv)1 SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0, 2 & pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out, 3 & Teta_out,u_out,ustar,tstar,L_mo,vhf,vvv) 4 4 IMPLICIT NONE 5 5 !======================================================================= … … 21 21 ! pz0(ngrid) surface roughness length 22 22 ! pg gravity (m s -2) 23 ! pz(ngrid,nlay) height of layers 23 ! zzlay(ngrid,nlay) height of mid-layers 24 ! zzlev(ngrid,nlay+1) height of layers interfaces 24 25 ! pu(ngrid,nlay) u component of the wind 25 26 ! pv(ngrid,nlay) v component of the wind 26 ! wmax(ngrid) maximum vertical velocity in thermals (might not be needed 27 ! if the computation of w* works) 27 ! wstar_in(ngrid) free convection velocity in thermals 28 28 ! hfmax(ngrid) maximum vertical turbulent heat flux in thermals 29 29 ! zmax(ngrid) height reached by the thermals (pbl height) … … 56 56 57 57 INTEGER, INTENT(IN) :: ngrid,nlay 58 REAL, INTENT(IN) :: pz0(ngrid) 59 REAL, INTENT(IN) :: pg, pz(ngrid,nlay)58 REAL, INTENT(IN) :: pz0(ngrid),ps(ngrid),pplay(ngrid,nlay) 59 REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay) 60 60 REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) 61 REAL, INTENT(IN) :: w max(ngrid),hfmax(ngrid),zmax(ngrid)61 REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) 62 62 REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) 63 63 REAL, INTENT(IN) :: z_out 64 65 ! Outputs: 66 ! -------- 67 64 68 REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid) 65 REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid),wstar(ngrid) 69 REAL T_out(ngrid) 70 REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) 71 REAL wstar(ngrid) 66 72 REAL, INTENT(OUT) :: L_mo(ngrid) 67 73 … … 139 145 ite=0. 140 146 residual=abs(pz0tcomp(ig)-pz0t) 141 zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1),wmax(ig)**2) 147 zu2(ig)=MAX(pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) 148 & ,(0.3*wstar_in(ig))**2) 142 149 143 150 DO WHILE((residual .gt. 0.01*pz0(ig)) .and. (ite .lt. 10.)) … … 147 154 ! Richardson number formulation proposed by D.E. England et al. (1995) 148 155 rib(ig) = (pg/ph(ig,1)) 149 & *sqrt( pz(ig,1)*pz0(ig))150 & *(((log( pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))156 & *sqrt(zzlev(ig,2)*pz0(ig)) 157 & *(((log(zzlev(ig,2)/pz0(ig)))**2)/(log(zzlev(ig,2)/pz0t))) 151 158 & *(ph(ig,1)-pts(ig))/zu2(ig) 152 159 ELSE … … 156 163 ENDIF 157 164 158 z1z0= pz(ig,1)/pz0(ig)159 z1z0t= pz(ig,1)/pz0t165 z1z0=zzlev(ig,2)/pz0(ig) 166 z1z0t=zzlev(ig,2)/pz0t 160 167 161 168 cdn(ig)=karman/log(z1z0) … … 310 317 ENDIF 311 318 319 T_out(:) = Teta_out(:)*(exp( 320 & (zout/zzlay(:,1))*(log(pplay(:,1)/ps)) 321 & ) 322 & )**rcp 323 324 312 325 !------------------------------------------------------------------------ 313 326 !------------------------------------------------------------------------ … … 325 338 DO k=1,nlay-1 326 339 DO ig=1, ngrid 327 IF (abs(zmax(ig)- pz(ig,k)) .lt.328 & abs(zmax(ig)-pz(ig,pbl_height_index(ig)))) THEN340 IF (abs(zmax(ig)-zzlay(ig,k)) .lt. 341 & abs(zmax(ig)-zzlay(ig,pbl_height_index(ig)))) THEN 329 342 pbl_height_index(ig)=k 330 343 ENDIF … … 337 350 DO k=1,nlay-1 338 351 DO ig=1, ngrid 339 dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(pz(ig,k+1)-pz(ig,k))352 dteta(ig,k) = (ph(ig,k+1)-ph(ig,k))/(zzlay(ig,k+1)-zzlay(ig,k)) 340 353 ENDDO 341 354 ENDDO … … 424 437 & 'Monin Obukhov length','m', 425 438 & dimout,L_mo) 426 call WRITEDIAGFI(ngrid,'w_star',427 & 'Free convection velocity','m',428 & dimout,wstar)429 call WRITEDIAGFI(ngrid,'z_i',430 & 'PBL height','m',431 & dimout,zmax)432 call WRITEDIAGFI(ngrid,'hf_max',433 & 'Maximum vertical heat flux','m',434 & dimout,hfmax)439 ! call WRITEDIAGFI(ngrid,'w_star', 440 ! & 'Free convection velocity','m', 441 ! & dimout,wstar) 442 ! call WRITEDIAGFI(ngrid,'z_i', 443 ! & 'PBL height','m', 444 ! & dimout,zmax) 445 ! call WRITEDIAGFI(ngrid,'hf_max', 446 ! & 'Maximum vertical heat flux','m', 447 ! & dimout,hfmax) 435 448 call WRITEDIAGFI(ngrid,'vvv', 436 449 & 'Vertical velocity variance at zout','m', -
trunk/LMDZ.MARS/libf/phymars/physiq.F
r528 r529 777 777 $ zdum1,zdum2,zdh,pdq,zflubid, 778 778 $ zdudif,zdvdif,zdhdif,zdtsdif,q2, 779 & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th) 779 & zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th 780 #ifdef MESOSCALE 781 & ,flag_LES 782 #endif 783 & ) 784 780 785 781 786 #ifdef MESOSCALE … … 1848 1853 if (calltherm .and. (z_out .gt. 0.)) then 1849 1854 1850 call pbl_parameters(ngrid,nlayer, z0,1851 & g,zzlay,z u,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,1855 call pbl_parameters(ngrid,nlayer,ps,zplay,z0, 1856 & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, 1852 1857 & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) 1853 1858 … … 1938 1943 ! THERMALS STUFF 1D 1939 1944 1940 z_out= 0.1945 z_out=1. 1941 1946 if (calltherm .and. (z_out .gt. 0.)) then 1942 1947 1943 call pbl_parameters(ngrid,nlayer, z0,1944 & g,zzlay,z u,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,1948 call pbl_parameters(ngrid,nlayer,ps,zplay,z0, 1949 & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, 1945 1950 & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) 1946 1951 -
trunk/LMDZ.MARS/libf/phymars/vdifc.F
r524 r529 5 5 $ pdufi,pdvfi,pdhfi,pdqfi,pfluxsrf, 6 6 $ pdudif,pdvdif,pdhdif,pdtsrf,pq2, 7 $ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax) 7 $ pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,hfmax 8 #ifdef MESOSCALE 9 & ,flag_LES 10 #endif 11 & ) 12 USE ioipsl_getincom 8 13 IMPLICIT NONE 9 14 … … 137 142 REAL qco2,mmean 138 143 144 #ifdef MESOSCALE 145 LOGICAL flag_LES !! pour LES avec isfflx!=0 146 #endif 139 147 c ** un petit test de coherence 140 148 c -------------------------- … … 242 250 DO ilev=1,nlay 243 251 DO ig=1,ngrid 244 qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep 252 qco2=MAX(1.E-30 253 & ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep) 245 254 c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) 246 255 mmean=1/(A*qco2 +B) … … 339 348 ENDIF 340 349 341 342 350 ! Some usefull diagnostics for the new surface layer parametrization : 343 344 ! 351 352 ! call WRITEDIAGFI(ngridmx,'pcdv', 345 353 ! & 'momentum drag','no units', 346 354 ! & 2,zcdv_true) … … 361 369 c ** schema de diffusion turbulente dans la couche limite 362 370 c ---------------------------------------------------- 363 ! 371 364 372 CALL vdif_kc(ptimestep,g,pzlev,pzlay 365 373 & ,pu,pv,ph,zcdv_true 366 374 & ,pq2,zkv,zkh,zq) 367 ! 375 368 376 ! pt(:,:)=ph(:,:)*ppopsk(:,:) 369 377 ! CALL yamada4(ngrid,nlay,ptimestep,g,r,pplev,pt … … 519 527 ENDDO 520 528 521 522 529 ! calcul de zc et zd pour la couche top en prenant en compte le terme 523 530 ! de variation de masse (on fait une boucle pour que ça converge) … … 526 533 527 534 llnt(:)=1 528 535 #ifdef MESOSCALE 536 IF (.not.flag_LES) THEN 537 #endif 529 538 DO ig=1,ngrid 530 539 DO l=1,nlay … … 537 546 5 continue 538 547 ENDDO 548 549 #ifdef MESOSCALE 550 ENDIF 551 #endif 539 552 540 553 DO ig=1,ngrid … … 589 602 ztsrf2(ig)=z1(ig)/z2(ig) 590 603 ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop 591 592 604 zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) 593 605 … … 633 645 c Using the wind modified by friction for lifting and sublimation 634 646 c ---------------------------------------------------------------- 635 c DO ig=1,ngrid 636 c zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) 637 c zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig)) 638 c zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig)) 639 c ENDDO 647 648 ! This is computed above and takes into account surface-atmosphere flux 649 ! enhancement by subgrid gustiness and atmospheric-stability related 650 ! variations of transfer coefficients. 651 652 ! DO ig=1,ngrid 653 ! zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) 654 ! zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig)) 655 ! zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig)) 656 ! ENDDO 640 657 641 658 c Calcul du flux vertical au bas de la premiere couche (dust) :
Note: See TracChangeset
for help on using the changeset viewer.
