AVIATION.f

Go to the documentation of this file.

 
!
      SUBROUTINE calllws(U,V,H,LLWS)
 
!
      USE vrbls2d, only: fis, u10, v10
      use params_mod, only: gi
      use ctlblk_mod, only: jsta, jend, im, jm, lsm, spval, ista, iend
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
!     DECLARE VARIABLES.
!     
      REAL,DIMENSION(IM,JM,LSM),INTENT(IN)    :: U,V,H
      REAL,DIMENSION(IM,JM),INTENT(INOUT)     :: LLWS
      REAL    :: Z1,Z2,HZ1,DH,U2,V2,W2,RT 
      INTEGER :: K1,K2 
      integer I,J,LP
 
!***************************************************************
!
!
 
      DO 100 j=jsta,jend
        DO i=ista,iend
 
          z1 = 10.0 + fis(i,j)*gi                              !Height of 10m levels geographic height (from sea level)
          
          IF(z1<h(i,j,lsm)) THEN                            !First search location of 10m wind level
            k1 = lsm + 1                                       !to see it is in which pressure levels
          ELSE
            DO lp = lsm,2,-1                                   !If not found, keep searching upward                              
             IF(z1>=h(i,j,lp).AND.z1<h(i,j,lp-1)) THEN
               k1 = lp 
             END IF
            END DO
          END IF
 
          hz1 = h(i,j,k1-1) - z1                                !Distance between K1-1 and 10m level
 
          dh = 0.0
 
          IF((hz1+10)>609.6) THEN                            !Then, search 2000ft(609.6m) location
            u2= u10(i,j) + (u(i,j,k1-1)-u10(i,j))*599.6/hz1     !found it between K1-1 and K1, then linear
            v2= v10(i,j) + (v(i,j,k1-1)-v10(i,j))*599.6/hz1     !interpolate to get wind at 2000ft U2,V2     
            z2= fis(i,j)*gi + 609.6
          ELSE                                                 !otherwise, keep on search upward
            DO lp = k1-1,2,-1
             dh=dh+(h(i,j,lp-1) - h(i,j,lp))
             IF((dh+hz1+10)>609.6) THEN                      !found the 2000ft level 
               z2=fis(i,j)*gi+609.6   
               rt=(z2-h(i,j,lp))/(h(i,j,lp-1)-h(i,j,lp))
               u2=u(i,j,lp)+rt*(u(i,j,lp-1)-u(i,j,lp))
               v2=v(i,j,lp)+rt*(v(i,j,lp-1)-v(i,j,lp))
               k2=lp
               exit
              END IF
             END DO
            END IF
 
!computer vector difference
         llws(i,j) = spval
         if(u10(i,j)<spval.and.v10(i,j)<spval)                   &
          llws(i,j)=sqrt((u2-u10(i,j))**2+(v2-v10(i,j))**2)/     &
                    609.6 * 1.943*609.6                         !unit: knot/2000ft
        ENDDO
 
100   CONTINUE     
 
      RETURN
      END
 
      SUBROUTINE calicing (T1,RH,OMGA, ICING)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .     
! SUBPROGRAM:    CALICING       COMPUTES In-Flight Icing
!   PRGRMMR: Binbin Zhou      /NCEP/EMC  DATE: 2005-08-16       
!     
! ABSTRACT:  
!    This program computes the in-flight icing condition
!    with the T-RH-OMGA algorithm provided by S. Silberberg of
!    NCEP/AWC (improved new version)
! 
!    According to S. Silberberg, Icing happens in following 
!    situation:
!       (1) -22C < T < 0C to      
!       (2)  RH > 70 %
!       (3) Ascent air, OMGA < 0 
!       (4) Equivalent Potential Vorticity (EPV) < 0
!       (5) Cloud water if SLD (supercooled large droplet)
!
!    Current version dosn't consider SLD, so cloud water           
!    is not used. EPV computation is not available for current
!    NCEP/EMC models(NAM, WRF, RSM), so EPV is also not
!    used
!
! USAGE:    CALL CALICING(T1,RH,OMGA,ICING)
!   INPUT ARGUMENT LIST:
!     T1     - TEMPERATURE (K)
!     RH     - RELATIVE HUMIDITY  (DECIMAL FORM)
!     OMGA   - Vertical velocity (Pa/sec)
!
!   OUTPUT ARGUMENT LIST: 
!     ICING     - ICING CONDITION (1 or 0)
!     
!   OUTPUT FILES:
!     NONE
!     
!   SUBPROGRAMS CALLED:
!     UTILITIES:
!     LIBRARY:
!       NONE
!     
!   ATTRIBUTES:
!     LANGUAGE: FORTRAN 90/77
!     MACHINE : BLUE AT NCEP
!$$$  
!
      use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!     
!     DECLARE VARIABLES.
!     
      REAL, DIMENSION(ista:iend,jsta:jend), INTENT(IN)    :: T1,RH,OMGA
      REAL, DIMENSION(ista:iend,jsta:jend), INTENT(INOUT) :: ICING 
      integer I,J
!***************************************************************
!
!
      DO j=jsta,jend
        DO i=ista,iend
        IF(omga(i,j)<spval.AND.t1(i,j)<spval.AND.rh(i,j)<spval) THEN
         IF(omga(i,j) < 0.0 .AND.                       &
            (t1(i,j) <= 273.0 .AND. t1(i,j) >= 251.0)   &
              .AND. rh(i,j) >= 70.0) THEN
 
           icing(i,j) = 1.0
         ELSE
           icing(i,j) = 0.0
         END IF
        ELSE
           icing(i,j) = spval
        ENDIF
        ENDDO
      ENDDO
 
      RETURN
      END
 
      SUBROUTINE calcat(U,V,H,U_OLD,V_OLD,H_OLD,CAT)
      use masks, only: dx, dy
      use ctlblk_mod, only: spval, jsta_2l, jend_2u, jsta_m, jend_m, &
              im, jm, ista_2l, iend_2u, ista_m, iend_m, ista, iend
      use gridspec_mod, only: gridtype
!
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
 
!     
!     DECLARE VARIABLES.
!     
      REAL,DIMENSION(ista_2l:iend_2u,jsta_2l:jend_2u),INTENT(IN)    :: U,V,H, &
                                                u_old,v_old,h_old
!      INTEGER,INTENT(IN)                      :: L
      REAL,DIMENSION(ista_2l:iend_2u,jsta_2l:jend_2u),INTENT(INOUT) :: CAT
 
      REAL  DSH, DST, DEF, CVG, VWS, TRBINDX
      INTEGER  IHE(JM),IHW(JM)
      integer I,J
      integer ISTART,ISTOP,JSTART,JSTOP
      real VWS1,VWS2,VWS3,VWS4
 
!***************************************************************
!
!
      cat=spval 
      DO j=jsta_2l,jend_2u
       IF(gridtype == 'A')THEN
        ihw(j)=-1
        ihe(j)=1 
    istart=ista_m
        istop=iend_m
        jstart=jsta_m
        jstop=jend_m
       ELSE IF(gridtype=='E')THEN
        ihw(j)=-mod(j,2)
        ihe(j)=ihw(j)+1
    istart=ista_m
        istop=iend_m
        jstart=jsta_m
        jstop=jend_m
       ELSE IF(gridtype=='B')THEN
        ihw(j)=-1
        ihe(j)=0 
    istart=ista_m
        istop=iend_m
        jstart=jsta_m
        jstop=jend_m
       ELSE 
        print*,'no gridtype specified, exit calcat comp'
    return  
       END IF   
      ENDDO
 
      call exch(u)
      call exch(v)
      call exch(u_old)
      call exch(v_old)
      call exch(h)
      call exch(h_old)
 
      DO 100 j=jstart,jstop
        DO i=istart,istop
!
          IF(gridtype=='B')THEN
           IF(u(i,j)<spval.and.u(i,j-1)<spval.and.u(i-1,j)<spval.and.u(i-1,j-1)<spval.and.&
              v(i,j)<spval.and.v(i,j-1)<spval.and.v(i-1,j)<spval.and.v(i-1,j-1)<spval)THEN
!dsh=dv/dx+du/dy 
           dsh=(0.5*(v(i,j)+v(i,j-1))-0.5*(v(i-1,j)+v(i-1,j-1)))*10000./dx(i,j) &
          +(0.5*(u(i,j)+u(i-1,j))-0.5*(u(i,j-1)+u(i-1,j-1)))*10000./dy(i,j)
!dst=du/dx-dv/dy
           dst =(0.5*(u(i,j)+u(i,j-1))-0.5*(u(i-1,j)+u(i-1,j-1)))*10000./dx(i,j) &
          -(0.5*(v(i,j)+v(i-1,j))-0.5*(v(i,j-1)+v(i-1,j-1)))*10000./dy(i,j)
           def = sqrt(dsh*dsh + dst*dst)
 
!cvg=-(du/dx+dv/dy)
           cvg = -((0.5*(u(i,j)+u(i,j-1))-0.5*(u(i-1,j)+u(i-1,j-1)))*10000./dx(i,j) &
                +(0.5*(v(i,j)+v(i-1,j))-0.5*(v(i,j-1)+v(i-1,j-1)))*10000./dy(i,j))     
           ELSE
            def = spval
            cvg = spval
           ENDIF
          ELSE
           IF(u(i,j+1)<spval.and.u(i,j-1)<spval.and.u(i+ihe(j),j)<spval.and.u(i+ihw(j),j)<spval.and.&
              v(i,j+1)<spval.and.v(i,j-1)<spval.and.v(i+ihe(j),j)<spval.and.v(i+ihw(j),j)<spval)THEN
!dsh=dv/dx+du/dy           
           dsh = (v(i+ihe(j),j) - v(i+ihw(j),j))*10000./(2*dx(i,j))   &  
              + (u(i,j+1) - u(i,j-1))*10000./(2*dy(i,j))
 
!dst=du/dx-dv/dy
           dst = (u(i+ihe(j),j) - u(i+ihw(j),j))*10000./(2*dx(i,j))   & 
              - (v(i,j+1) - v(i,j-1))*10000./(2*dy(i,j))
 
           def = sqrt(dsh*dsh + dst*dst)
 
!cvg=-(du/dx+dv/dy)
           cvg = -( (u(i+ihe(j),j) - u(i+ihw(j),j))*10000./(2*dx(i,j)) &
                  +(v(i,j+1) - v(i,j-1))*10000./(2*dy(i,j)) )
           ELSE
            def = spval
            cvg = spval
           ENDIF
          END IF
      
          IF(gridtype == 'A')THEN
!vws=d|U|/dz
           IF(u_old(i,j)<spval.and.u(i,j)<spval.and.&
              v_old(i,j)<spval.and.v(i,j)<spval.and.&
              h_old(i,j)<spval.and.h(i,j)<spval)THEN
           vws = ( sqrt(u_old(i,j)**2+v_old(i,j)**2 ) -               &
                  sqrt(u(i,j)**2+v(i,j)**2 )   ) *                    &
                  1000.0/(h_old(i,j) - h(i,j))
           ELSE
            vws = spval
           ENDIF
          else IF(gridtype == 'E')THEN
!vws=d|U|/dz
           IF(u_old(i+ihe(j),j)<spval.and.u(i+ihe(j),j)<spval.and.&
              v_old(i+ihe(j),j)<spval.and.v(i+ihe(j),j)<spval)THEN
            
       vws1 = ( sqrt(u_old(i+ihe(j),j)**2+v_old(i+ihe(j),j)**2 ) -&
                  sqrt(u(i+ihe(j),j)**2+v(i+ihe(j),j)**2 )   ) 
           ELSE
            vws1 = spval
           ENDIF
!vws=d|U|/dz
           IF(u_old(i+ihw(j),j)<spval.and.u(i+ihw(j),j)<spval.and.&
              v_old(i+ihw(j),j)<spval.and.v(i+ihw(j),j)<spval)THEN
           vws2 = ( sqrt(u_old(i+ihw(j),j)**2+v_old(i+ihw(j),j)**2 ) -&   
                  sqrt(u(i+ihw(j),j)**2+v(i+ihw(j),j)**2 )   ) 
           ELSE
            vws2 = spval
           ENDIF
!vws=d|U|/dz
           IF(u_old(i,j-1)<spval.and.u(i,j-1)<spval.and.&
              v_old(i,j-1)<spval.and.v(i,j-1)<spval)THEN
           vws3 = ( sqrt(u_old(i,j-1)**2+v_old(i,j-1)**2 ) -          & 
                  sqrt(u(i,j-1)**2+v(i,j-1)**2 )   ) 
           ELSE
            vws3 = spval
           ENDIF
!vws=d|U|/dz
           IF(u_old(i,j+1)<spval.and.u(i,j+1)<spval.and.&
              v_old(i,j+1)<spval.and.v(i,j+1)<spval)THEN
           vws4 = ( sqrt(u_old(i,j+1)**2+v_old(i,j+1)**2 ) -          & 
                  sqrt(u(i,j+1)**2+v(i,j+1)**2 )   ) 
           ELSE
            vws4 = spval
           ENDIF
 
           IF(vws1<spval.and.vws2<spval.and.vws3<spval.and.vws4<spval.and.&
              h_old(i,j)<spval.and.h(i,j)<spval)THEN
           vws=1000.0*(vws1+vws2+vws3+vws4)/4.0/(h_old(i,j) - h(i,j))
           ELSE
            vws = spval
           ENDIF
      ELSE IF(gridtype == 'B')THEN
           IF(u_old(i+ihe(j),j)<spval.and.u(i+ihe(j),j)<spval.and.&
              v_old(i+ihe(j),j)<spval.and.v(i+ihe(j),j)<spval)THEN
       vws1 = ( sqrt(u_old(i+ihe(j),j)**2+v_old(i+ihe(j),j)**2 ) -&
                  sqrt(u(i+ihe(j),j)**2+v(i+ihe(j),j)**2 )   ) 
           ELSE
            vws1 = spval
           ENDIF
!vws=d|U|/dz
           IF(u_old(i+ihw(j),j)<spval.and.u(i+ihw(j),j)<spval.and.&
              v_old(i+ihw(j),j)<spval.and.v(i+ihw(j),j)<spval)THEN
           vws2 = ( sqrt(u_old(i+ihw(j),j)**2+v_old(i+ihw(j),j)**2 ) -&   
                  sqrt(u(i+ihw(j),j)**2+v(i+ihw(j),j)**2 )   ) 
           ELSE
            vws2 = spval
           ENDIF
!vws=d|U|/dz
           IF(u_old(i,j-1)<spval.and.u(i,j-1)<spval.and.&
              v_old(i,j-1)<spval.and.v(i,j-1)<spval)THEN
           vws3 = ( sqrt(u_old(i,j-1)**2+v_old(i,j-1)**2 ) -          & 
                  sqrt(u(i,j-1)**2+v(i,j-1)**2 )   ) 
           ELSE
            vws3 = spval
           ENDIF
!vws=d|U|/dz
           IF(u_old(i-1,j-1)<spval.and.u(i-1,j-1)<spval.and.&
              v_old(i-1,j-1)<spval.and.v(i-1,j-1)<spval)THEN
           vws4 = ( sqrt(u_old(i-1,j-1)**2+v_old(i-1,j-1)**2 ) -          & 
                  sqrt(u(i-1,j-1)**2+v(i-1,j-1)**2 )   ) 
           ELSE
            vws4 = spval
           ENDIF
 
           IF(vws1<spval.and.vws2<spval.and.vws3<spval.and.vws4<spval.and.&
              h_old(i,j)<spval.and.h(i,j)<spval)THEN
           vws=1000.0*(vws1+vws2+vws3+vws4)/4.0/(h_old(i,j) - h(i,j)) 
           ELSE
            vws=spval
           ENDIF
      END IF  
          
         IF(vws<spval.and.def<spval.and.cvg<spval)THEN 
          trbindx = abs(vws)*(def + abs(cvg))
      
          IF(trbindx<=4.) THEN
            cat(i,j) = 0.0
          ELSE IF(trbindx<=8.) THEN
            cat(i,j)=1.0
          ELSE IF(trbindx<=12.) THEN
            cat(i,j)=2.0
          ELSE
            cat(i,j)=3.0
          END IF        
         ELSE
          cat(i,j)=spval
         ENDIF
        ENDDO
 
100   CONTINUE     
 
      RETURN
      END
 
      SUBROUTINE calceiling (CLDZ,TCLD,CEILING)
      USE vrbls2d, only: fis
      use params_mod, only: small, gi
      use ctlblk_mod, only: jsta, jend, spval, im, modelname, ista, iend
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!     
!     DECLARE VARIABLES.
!     
      REAL, DIMENSION(ista:iend,jsta:jend), INTENT(IN)    :: CLDZ, TCLD
      REAL, DIMENSION(ista:iend,jsta:jend), INTENT(INOUT) :: CEILING
      integer I,J
!***************************************************************
!
!
      DO j=jsta,jend
        DO i=ista,iend
          IF(abs(tcld(i,j)-spval) <= small) THEN
            ceiling(i,j)=spval
          ELSE IF(tcld(i,j) >= 50.) THEN
            if(modelname == 'RAPR')then
              ceiling(i,j) = cldz(i,j) - fis(i,j)*gi
            else
              ceiling(i,j) = cldz(i,j) ! for RAP/HRRR   - FIS(I,J)*GI
            endif
          ELSE
            ceiling(i,j) = 20000.0
          END IF
 
          IF(ceiling(i,j) < 0.0) ceiling(i,j)=20000.0
 
        ENDDO
      ENDDO
 
      RETURN
      END
 
      SUBROUTINE calfltcnd (CEILING,FLTCND)
      use vrbls2d, only: vis
      use ctlblk_mod, only: jsta, jend, im, spval, ista, iend
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!     
!     DECLARE VARIABLES.
!     
      REAL, DIMENSION(ista:iend,jsta:jend), INTENT(IN)    :: CEILING
      REAL, DIMENSION(ista:iend,jsta:jend), INTENT(INOUT) :: FLTCND
      REAL  CEIL,VISI
      integer I,J
!
!***************************************************************
!
!
      DO j=jsta,jend
        DO i=ista,iend
 
        IF(ceiling(i,j)<spval.and.vis(i,j)<spval)THEN
          ceil = ceiling(i,j) * 3.2808               !from m -> feet
          visi = vis(i,j) / 1609.0                   !from m -> miles       
 
          IF(ceil<500.0 .OR. visi<1.0 ) THEN
             fltcnd(i,j) = 1.0
 
          ELSE IF( (ceil>=500.AND.ceil<1000.0) .OR.          &
                   (visi>=1.0.AND.visi<3.0) ) THEN
             fltcnd(i,j) = 2.0
 
          ELSE IF( (ceil>=1000.AND.ceil<=3000.0) .OR.         &
                   (visi>=3.0.AND.visi<=5.0) ) THEN
             fltcnd(i,j) = 3.0
 
          ELSE IF( ceil>3000.0  .OR. visi>5.0) THEN
             fltcnd(i,j) = 4.0
 
          END IF
        ELSE
          fltcnd(i,j) = spval
        ENDIF
        ENDDO
      ENDDO
 
      RETURN
      END