CALGUST.f

Go to the documentation of this file.

!

      SUBROUTINE calgust(LPBL,ZPBL,GUST)

!     
!     
      use vrbls3d, only: uh, vh, zint, zmid
      use vrbls2d , only: u10h, v10h, u10,v10, fis
      use params_mod, only: d25, gi
      use ctlblk_mod, only: jsta, jend, spval, jsta_m, jend_m, num_procs, mpi_comm_comp, lm,&
              modelname, im, jm, jsta_2l, jend_2u, ista, iend, ista_m, iend_m, ista_2l, iend_2u
      use gridspec_mod, only: gridtype

      implicit none

      include "mpif.h"
! 
!     INCLUDE ETA GRID DIMENSIONS.  SET/DERIVE PARAMETERS.
!
!     DECLARE VARIABLES.
!     
      INTEGER,intent(in) :: lpbl(ista_2l:iend_2u,jsta_2l:jend_2u)
      REAL,intent(in)    :: zpbl(ista_2l:iend_2u,jsta_2l:jend_2u)
      REAL,intent(inout) :: gust(ista_2l:iend_2u,jsta_2l:jend_2u)

      integer i,j,ie,iw, l, k, istart, istop, jstart, jstop
      integer lmin,lxxx,ierr
      real zsfc,delwind,usfc,vsfc,sfcwind,wind,u0,v0,dz
!     
!     
!*****************************************************************************
!     START CALMXW HERE.
!     
!     LOOP OVER THE GRID.
!    
!$omp parallel do private(i,j)
      DO j=jsta,jend
        DO i=ista,iend
          gust(i,j) = spval 
        ENDDO
      ENDDO
      
      IF(gridtype == 'A') THEN
        istart = ista
        istop  = iend
        jstart = jsta
        jstop  = jend
      ELSE
        istart = ista_m
        istop  = iend_m
        jstart = jsta_m
        jstop  = jend_m
        if ( num_procs > 1 ) then
         !CALL EXCH(U10(1,jsta_2l))
         !CALL EXCH(V10(1,jsta_2l))
         lmin = max(1, minval(lpbl(ista:iend,jsta:jend)))
         CALL mpi_allreduce(lmin,lxxx,1,mpi_integer,mpi_min,mpi_comm_comp,ierr)
         DO l=lxxx,lm
          CALL exch(uh(1,jsta_2l,l))
          CALL exch(vh(1,jsta_2l,l))
         END DO 
        END IF 
      END IF
!
!     ASSUME THAT U AND V HAVE UPDATED HALOS
!
!!$omp  parallel do private(i,j,ie,iw,mxww,u0,v0,wind)
      DO j=jstart,jstop
        DO i=istart,istop
           l=lpbl(i,j) 
           IF(gridtype == 'E') THEN
             ie = i + mod(j+1,2) 
             iw = i + mod(j+1,2)-1

           if(u10h(i,j)<spval.and.uh(i,j+1,l)<spval.and.uh(ie,j,l)<spval.and.uh(iw,j,l)<spval.and.uh(i,j-1,l)<spval) then

!        USFC=D25*(U10(I,J-1)+U10(IW,J)+U10(IE,J)+U10(I,J+1)) 
!        VSFC=D25*(V10(I,J-1)+V10(IW,J)+V10(IE,J)+V10(I,J+1))
             usfc = u10h(i,j)
             vsfc = v10h(i,j)
             sfcwind = sqrt(usfc*usfc + vsfc*vsfc)
             u0   = d25*(uh(i,j-1,l)+uh(iw,j,l)+uh(ie,j,l)+uh(i,j+1,l))
             v0   = d25*(vh(i,j-1,l)+vh(iw,j,l)+vh(ie,j,l)+vh(i,j+1,l))
             wind = sqrt(u0*u0 + v0*v0)

             else
             wind = spval
             endif

           ELSE IF(gridtype == 'B') THEN
             ie = i 
             iw = i-1

!        USFC=D25*(U10(I,J-1)+U10(IW,J)+U10(IE,J)+U10(IW,J-1)) 
!        VSFC=D25*(V10(I,J-1)+V10(IW,J)+V10(IE,J)+V10(IW,J-1))

             if(u10h(i,j)<spval.and.uh(iw,j-1,l)<spval) then

             usfc = u10h(i,j)
             vsfc = v10h(i,j)
             sfcwind = sqrt(usfc*usfc + vsfc*vsfc)
             u0 = d25*(uh(i,j-1,l)+uh(iw,j,l)+uh(ie,j,l)+uh(iw,j-1,l))
             v0 = d25*(vh(i,j-1,l)+vh(iw,j,l)+vh(ie,j,l)+vh(iw,j-1,l))
             wind = sqrt(u0*u0 + v0*v0) 
             else
             wind = spval
             endif
           ELSE IF(gridtype == 'A') THEN

             usfc    = u10(i,j)
             vsfc    = v10(i,j)
             if (usfc < spval .and. vsfc < spval) then
               sfcwind = sqrt(usfc*usfc + vsfc*vsfc) 
             else
               sfcwind = spval
             endif
             if(modelname == 'RAPR') then
               zsfc = zint(i,j,lm+1)
               l = lpbl(i,j)
! in RUC do 342 k=2,k1-1, where k1 - first level above PBLH
               gust(i,j) = sfcwind
               do k=lm-1,l-1,-1

                 if(uh(i,j,l)<spval) then
                 u0   = uh(i,j,k)
                 v0   = vh(i,j,k)
                 wind = sqrt(u0*u0 + v0*v0)
                 delwind = wind - sfcwind
                 dz = zmid(i,j,k)-zsfc
                 delwind  = delwind*(1.0-min(0.5,dz/2000.))
                 gust(i,j) = max(gust(i,j),sfcwind+delwind)
               else
                 gust(i,j) = spval
               endif
               enddo
             else
               if(uh(i,j,l)<spval) then
               u0   = uh(i,j,l)
               v0   = vh(i,j,l)
               wind = sqrt(u0*u0 + v0*v0 )
               else
               wind = spval
               endif
             endif ! endif RAPR

           ELSE
!            print*,'unknown grid type, not computing wind gust'
             return
           END IF

           if(modelname /= 'RAPR')then
             if (sfcwind < spval) then
               delwind   = wind - sfcwind
               zsfc      = fis(i,j)*gi
               delwind   = delwind*(1.0-min(0.5,zpbl(i,j)/2000.))
               gust(i,j) = sfcwind + delwind
             else
               gust(i,j) = wind
             endif
           endif
        enddo
      enddo

!     END OF ROUTINE.
!     
      RETURN
      END