Ignore:
Timestamp:
Apr 27, 2012, 11:22:46 AM (13 years ago)
Author:
acolaitis
Message:

Updated Monin-Obukhov interpolation scheme: corrected early afternoon interpolations problem at peak convective activity, corrected nightime interpolation problems in tricky conditions around the critical richardson number.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F

    r605 r636  
    144144      lambda=(sqrt(bh/bm))/alphah
    145145      ric=betah/(betam**2)
    146 
    147146      DO ig=1,ngrid
    148147       ite=0.
    149148       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
    152153
    153154       DO WHILE((residual .gt. 0.01*pz0(ig)) .and.  (ite .lt. 10.))
    … …  
    157158            ! Richardson number formulation proposed by D.E. England et al. (1995)
    158159          rib(ig) = (pg/pts(ig))
    159      &        *sqrt(zzlev(ig,2)*pz0(ig))
    160      &        *(((log(zzlev(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)))
    161162     &        *(ph(ig,1)-pts(ig))/zu2(ig)
    162163         ELSE
    … …  
    166167         ENDIF
    167168
    168          z1z0=zzlev(ig,2)/pz0(ig)
    169          z1z0t=zzlev(ig,2)/pz0t
     169         z1z0=zzlay(ig,1)/pz0(ig)
     170         z1z0t=zzlay(ig,1)/pz0t
    170171
    171172         cdn(ig)=karman/log(z1z0)
    … …  
    221222
    222223
    223 ! Large-scale wind at first layer (without gustiness) and
    224224! u* theta* computation
     225! and Monin Obukhov length:
    225226
    226227      DO ig=1,ngrid
    … …  
    228229           ustar(ig)=0.
    229230           tstar(ig)=0.
     231           L_mo(ig)=0.
    230232         ELSE
    231233           ustar(ig)=sqrt(pcdv(ig))
    … …  
    233235           tstar(ig)=-pcdh(ig)*(pts(ig)-ph(ig,1))
    234236     &        /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
    235238         ENDIF
    236239      ENDDO
    … …  
    238241! Monin Obukhov length:
    239242
    240       DO ig=1,ngrid
    241          IF (rib(ig) .gt. ric) THEN
    242             L_mo(ig)=0.
    243          ELSE
    244             L_mo(ig)=pts(ig)*(ustar(ig)**2)/(pg*karman*tstar(ig))  !as defined here, L is positive for SBL, negative for UBL
    245          ENDIF
    246       ENDDO
     243!      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
    247250
    248251! Interpolation:
    … …  
    253256           Teta_out(ig,n)=pts(ig)
    254257        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
    261307        ELSEIF (L_mo(ig) .lt. 0.) THEN
    262308
    … …  
    304350          ENDIF
    305351        ELSE
    306            u_out(ig,n)=0.
    307            Teta_out(ig,n)=pts(ig)
     352           u_out(ig,n)=sqrt(pu(ig,1)**2 + pv(ig,1)**2)
     353           Teta_out(ig,n)=ph(ig,1)
    308354        ENDIF
    309355        IF(zout .lt. pz0(ig)) THEN
    310356           u_out(ig,n)=0.
    311357        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
    312378      ENDDO
    313379
    … …  
    420486      ENDIF ! of if calltherm
    421487
     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
    422508      RETURN
    423509      END
Note: See TracChangeset for help on using the changeset viewer.