CALTAU.f

Go to the documentation of this file.

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
      SUBROUTINE caltau(TAUX,TAUY)

!     
!
      use vrbls3d,      only: zint, pmid, q, t, uh, vh, el_pbl, zmid
      use vrbls2d,      only: z0, uz0, vz0
      use masks,        only: lmh
      use params_mod,   only:  d00, d50, h1, d608, rd, d25
      use ctlblk_mod,   only: jsta_2l, jend_2u, lm, jsta, jend, spval, jsta_m,&
                            jm, im, jend_m, ista, iend, ista_m, iend_m, ista_2l, iend_2u
      use gridspec_mod, only: gridtype

      implicit none
!
!     DECLARE VARIABLES.
      INTEGER, dimension(4)  :: kk(4)
      INTEGER, dimension(jm) :: ive, ivw
      REAL, dimension(ista:iend,jsta:jend), intent(inout) :: taux, tauy
      REAL, ALLOCATABLE :: el(:,:,:)
      REAL, dimension(ista:iend,jsta:jend) ::  egridu,egridv,egrid4,egrid5, el0
      REAL uz0v,vz0v
      CHARACTER*1 agrid
      integer i,j,lmhk,ie,iw,ii,jj
      real dz,rdz,rsfc,tv,rho,ulmh,vlmh,deludz,delvdz,elsqr,zint1,    &
           zint2,z0v,psfc,tvv,qvv,elv,elv1,elv2
!     
!********************************************************************
!     START CALTAU HERE.
!    
      ALLOCATE (el(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
!
!     COMPUTE MASTER LENGTH SCALE.
!
!      CALL CLMAX(EL0,EGRIDU,EGRIDV,EGRID4,EGRID5)
!      CALL MIXLEN(EL0,EL)
!     
!     INITIALIZE OUTPUT AND WORK ARRAY TO ZERO.
!     
      DO j=jsta,jend
      DO i=ista,iend
        egridu(i,j) = d00
        egridv(i,j) = d00
        taux(i,j)   = spval
        tauy(i,j)   = spval
      ENDDO
      ENDDO
!     
!     COMPUTE SURFACE LAYER U AND V WIND STRESSES.
!
!     ASSUME THAT U AND V HAVE UPDATED HALOS
!
      IF(gridtype == 'A')THEN
       CALL clmax(el0,egridu,egridv,egrid4,egrid5)
       CALL mixlen(el0,el)

       DO j=jsta,jend
       DO i=ista,iend
!
        lmhk = nint(lmh(i,j))
        IF(el(i,j,lmhk-1)<spval.and.z0(i,j)<spval.and. &
           uz0(i,j)<spval.and.vz0(i,j)<spval)THEN
!
!       COMPUTE THICKNESS OF LAYER AT MASS POINT.
!
        dz  = d50*(zint(i,j,lmhk)-zint(i,j,lmhk+1))
        dz  = dz-z0(i,j)
        rdz = 1./dz
!
!        COMPUTE REPRESENTATIVE AIR DENSITY.
!
        psfc = pmid(i,j,lmhk)
        tv   = (h1+d608*q(i,j,lmhk))*t(i,j,lmhk)
        rho  = psfc/(rd*tv)
!     
!        COMPUTE A MEAN MASS POINT WIND IN THE 
!        FIRST ATMOSPHERIC ETA LAYER.
!
        ulmh = uh(i,j,lmhk)
        vlmh = vh(i,j,lmhk)
!
!       COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
!
        deludz = (ulmh-uz0(i,j))*rdz
        delvdz = (vlmh-vz0(i,j))*rdz
!     
!       COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
!
        elsqr     = el(i,j,lmhk-1)*el(i,j,lmhk-1)
        taux(i,j) = rho*elsqr*deludz*deludz
        tauy(i,j) = rho*elsqr*delvdz*delvdz
        ELSE
        taux(i,j) = spval
        tauy(i,j) = spval
        ENDIF

!
       END DO
       END DO
      ELSE IF(gridtype == 'E')THEN
       call exch(zint(1,jsta_2l,lm))
       call exch(zint(1,jsta_2l,lm+1))
       call exch(z0(1,jsta_2l))
       call exch(pmid(1,jsta_2l,lm))
       call exch(t(1,jsta_2l,lm))
       call exch(q(1,jsta_2l,lm))
       call exch(el_pbl(1,jsta_2l,lm))
       call exch(el_pbl(1,jsta_2l,lm-1))

       DO j=jsta_m,jend_m
        ive(j)=mod(j,2)
        ivw(j)=ive(j)-1
       ENDDO
 
       DO j=jsta_m,jend_m
       DO i=ista_m,iend_m
!
        lmhk = nint(lmh(i,j)) 
        ie=i+ive(j)
        iw=i+ivw(j)
        zint1=(zint(iw,j,lmhk)+zint(ie,j,lmhk)                  &  
         +zint(i,j+1,lmhk)+zint(i,j-1,lmhk))*d25
        zint2=(zint(iw,j,lmhk+1)+zint(ie,j,lmhk+1)              &
         +zint(i,j+1,lmhk+1)+zint(i,j-1,lmhk+1))*d25
        dz  = d50*(zint1-zint2)       
        z0v=(z0(iw,j)+z0(ie,j)+z0(i,j+1)+z0(i,j-1))*d25
        dz  = dz-z0v
        rdz = 1./dz
!
!        COMPUTE REPRESENTATIVE AIR DENSITY.
!
        psfc = (pmid(iw,j,lmhk)+pmid(ie,j,lmhk)                &
         +pmid(i,j+1,lmhk)+pmid(i,j-1,lmhk))*d25 
        tvv = (t(iw,j,lmhk)+t(ie,j,lmhk)                       &
         +t(i,j+1,lmhk)+t(i,j-1,lmhk))*d25
        qvv = (q(iw,j,lmhk)+q(ie,j,lmhk)                       &
         +q(i,j+1,lmhk)+q(i,j-1,lmhk))*d25
        tv   = (h1+d608*qvv)*tvv
        rho  = psfc/(rd*tv) 

!       COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
!
        deludz = (uh(i,j,lmhk)-uz0(i,j))*rdz
        delvdz = (vh(i,j,lmhk)-vz0(i,j))*rdz 

!       COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
!                                       
        elv1=(el_pbl(iw,j,lmhk)+el_pbl(ie,j,lmhk)              &
         +el_pbl(i,j+1,lmhk)+el_pbl(i,j-1,lmhk))*d25
        elv2=(el_pbl(iw,j,lmhk-1)+el_pbl(ie,j,lmhk-1)          &
         +el_pbl(i,j+1,lmhk-1)+el_pbl(i,j-1,lmhk-1))*d25
        elv=(elv1+elv2)/2.0  ! EL is defined at the bottom of layer
        elsqr       =elv*elv
        taux(i,j)=rho*elsqr*deludz*deludz 
        tauy(i,j)=rho*elsqr*delvdz*delvdz
!   ii=im/2
!   jj=(jsta+jend)/2
!        if(i==ii.and.j==jj)print*,'sample tau'
!     & ,RHO,ELSQR,DELUDZ,DELVDZ  
       END DO
       END DO
      ELSE IF(gridtype == 'B')THEN
! PUT TAUX AND TAUY ON MASS POINTS      
       call exch(vh(1,jsta_2l,lm))
       DO j=jsta_m,jend_m
       DO i=ista_m,iend_m
!
        lmhk = nint(lmh(i,j))
!
!       COMPUTE THICKNESS OF LAYER AT MASS POINT.
!
!        DZ  = D50*(ZINT(I,J,LMHK)-ZINT(I,J,LMHK+1))
!        DZ  = ZMID(I,J,LMHK)-Z0(I,J)
        dz=zmid(i,j,lmhk)-(z0(i,j)+zint(i,j,lmhk+1))
        if(dz==0.0)dz=0.2
        rdz = 1./dz
!
!        COMPUTE REPRESENTATIVE AIR DENSITY.
!
        psfc = pmid(i,j,lmhk)
        tv   = (h1+d608*q(i,j,lmhk))*t(i,j,lmhk)
        rho  = psfc/(rd*tv)
!     
!        PUT U AND V ONTO MASS POINTS
!
        ulmh = 0.5*(uh(i-1,j,lmhk)+uh(i,j,lmhk))
        vlmh = 0.5*(vh(i,j-1,lmhk)+vh(i,j,lmhk))
!
!       COMPUTE WIND SHEAR COMPONENTS ACROSS LAYER.
!
        deludz = (ulmh-uz0(i,j))*rdz
        delvdz = (vlmh-vz0(i,j))*rdz
!     
!       COMPUTE U (EGRIDU) AND V (EGRIDV) WIND STRESSES.
!
        elv=0.5*(el_pbl(i,j,lmhk)+el_pbl(i,j,lmhk-1))
        elsqr     = elv*elv
        taux(i,j) = rho*elsqr*deludz*deludz
!        if(TAUX(I,J)>1.0e2)print*,'Debug TAUX= ',i,j, &
!       ELV,ULMH,UZ0(I,J),ZMID(I,J,LMHK),Z0(I,J),RDZ,TAUX(I,J),zint(i,j,lm+1)
        tauy(i,j) = rho*elsqr*delvdz*delvdz
    
       END DO
       END DO
      END IF
!     
      DEALLOCATE(el)
!     END OF ROUTINE.
!     
      RETURN
      END