velo.f90

FDS为火灾动力学模拟软件源代码,该软件为开源项目,代码语言主要为FORTRAN,可在WINDOWS和LINUX下编译运行,详细说明可参考http://fire.nist.gov/fds/官方网址

   IIO2 = IJKE(12,IE)
   JJO2 = IJKE(13,IE)
   KKO2 = IJKE(14,IE)
   SF  => SURFACE(IBC)
   BC  = SF%SLIP_FACTOR
   IF (BC>1.5_EB) THEN
      IF (SF%T_IGN==T_BEGIN) THEN
         TSI = T
      ELSE
         TSI=T-SF%T_IGN
      ENDIF
      FVT = EVALUATE_RAMP(TSI,SF%TAU(TIME_VELO),SF%RAMP_INDEX(TIME_VELO))
   ENDIF
 
   COMPONENT: SELECT CASE(IEC)
 
      CASE(1) COMPONENT

         ICMM = CELL_INDEX(II,JJ,KK)
         ICMP = CELL_INDEX(II,JJ,KK+1)
         ICPM = CELL_INDEX(II,JJ+1,KK)
         ICPP = CELL_INDEX(II,JJ+1,KK+1)
         VP   = VV(II,JJ,KK+1)
         VM   = VV(II,JJ,KK)
         WP   = WW(II,JJ+1,KK)
         WM   = WW(II,JJ,KK)

         IF (NOM1==NM) THEN
            IWM = WALL_INDEX(ICMM, 3) 
            IWP = WALL_INDEX(ICPM, 3) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) VP = BC*VM
               IF (BC>1.5_EB) VP = FVT*SF%VEL_T(2)
            ENDIF
         ENDIF
         IF (NOM1==-NM) THEN
            IWM = WALL_INDEX(ICMP,-3) 
            IWP = WALL_INDEX(ICPP,-3) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) VM = BC*VP
               IF (BC>1.5_EB) VM = FVT*SF%VEL_T(2)
            ENDIF
         ENDIF

         IF (NOM2==NM) THEN
            IWM = WALL_INDEX(ICMM, 2) 
            IWP = WALL_INDEX(ICMP, 2) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) WP = BC*WM
               IF (BC>1.5_EB) WP = FVT*SF%VEL_T(2)
            ENDIF
         ENDIF
         IF (NOM2==-NM) THEN
            IWM = WALL_INDEX(ICPM,-2) 
            IWP = WALL_INDEX(ICPP,-2) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) WM = BC*WP
               IF (BC>1.5_EB) WM = FVT*SF%VEL_T(2)
            ENDIF
         ENDIF

         IF (ABS(NOM1)/=NM .AND. NOM1/=0) THEN
            IF (NOM1<0) VM = OMESH(ABS(NOM1))%V(IIO1,JJO1,KKO1)
            IF (NOM1>0) VP = OMESH(ABS(NOM1))%V(IIO1,JJO1,KKO1)
         ENDIF
         IF (ABS(NOM2)/=NM .AND. NOM2/=0) THEN
            IF (NOM2<0) WM = OMESH(ABS(NOM2))%W(IIO2,JJO2,KKO2)
            IF (NOM2>0) WP = OMESH(ABS(NOM2))%W(IIO2,JJO2,KKO2)
         ENDIF

         MUA = .25_EB*( MU(II,JJ,KK) + MU(II,JJ+1,KK) + MU(II,JJ+1,KK+1) + MU(II,JJ,KK+1) )
         DVDZ = RDZN(KK)*(VP-VM)
         DWDY = RDYN(JJ)*(WP-WM)
         OME_E(IE) = DWDY - DVDZ
         TAU_E(IE) = MUA*(DVDZ + DWDY)
         IF (JJ==0)    WW(II,JJ,KK)   = WM
         IF (JJ==JBAR) WW(II,JJ+1,KK) = WP
         IF (KK==0)    VV(II,JJ,KK)   = VM
         IF (KK==KBAR) VV(II,JJ,KK+1) = VP
         IF (CORRECTOR .AND. JJ>0 .AND. JJ<JBAR .AND. KK>0 .AND. KK<KBAR) THEN
            W_Y(II,JJ,KK) = 0.5_EB*(WM+WP)
            V_Z(II,JJ,KK) = 0.5_EB*(VM+VP)
         ENDIF
 
      CASE(2) COMPONENT

         ICMM = CELL_INDEX(II,JJ,KK)
         ICMP = CELL_INDEX(II,JJ,KK+1)
         ICPM = CELL_INDEX(II+1,JJ,KK)
         ICPP = CELL_INDEX(II+1,JJ,KK+1)
         UP   = UU(II,JJ,KK+1)
         UM   = UU(II,JJ,KK)
         WP   = WW(II+1,JJ,KK)
         WM   = WW(II,JJ,KK)

         IF (NOM1==NM) THEN
            IWM = WALL_INDEX(ICMM, 3) 
            IWP = WALL_INDEX(ICPM, 3) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) UP = BC*UM
               IF (BC>1.5_EB) UP = FVT*SF%VEL_T(1)
            ENDIF
         ENDIF
         IF (NOM1==-NM) THEN
            IWM = WALL_INDEX(ICMP,-3) 
            IWP = WALL_INDEX(ICPP,-3) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) UM = BC*UP
               IF (BC>1.5_EB) UM = FVT*SF%VEL_T(1)
            ENDIF
         ENDIF

         IF (NOM2==NM) THEN
            IWM = WALL_INDEX(ICMM, 1) 
            IWP = WALL_INDEX(ICMP, 1) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) WP = BC*WM
               IF (BC>1.5_EB) WP = FVT*SF%VEL_T(2)
            ENDIF
         ENDIF
         IF (NOM2==-NM) THEN
            IWM = WALL_INDEX(ICPM,-1) 
            IWP = WALL_INDEX(ICPP,-1) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) WM = BC*WP
               IF (BC>1.5_EB) WM = FVT*SF%VEL_T(2)
            ENDIF
         ENDIF

         IF (ABS(NOM1)/=NM .AND. NOM1/=0) THEN
            IF (NOM1<0) UM = OMESH(ABS(NOM1))%U(IIO1,JJO1,KKO1)
            IF (NOM1>0) UP = OMESH(ABS(NOM1))%U(IIO1,JJO1,KKO1)
         ENDIF
         IF (ABS(NOM2)/=NM .AND. NOM2/=0) THEN
            IF (NOM2<0) WM = OMESH(ABS(NOM2))%W(IIO2,JJO2,KKO2)
            IF (NOM2>0) WP = OMESH(ABS(NOM2))%W(IIO2,JJO2,KKO2)
         ENDIF

         MUA = .25_EB*( MU(II,JJ,KK) + MU(II+1,JJ,KK) + MU(II+1,JJ,KK+1) + MU(II,JJ,KK+1) )
         DUDZ = RDZN(KK)*(UP-UM)
         DWDX = RDXN(II)*(WP-WM)
         OME_E(IE) = DUDZ - DWDX
         TAU_E(IE) = MUA*(DUDZ + DWDX)
         IF (II==0)    WW(II,JJ,KK)   = WM
         IF (II==IBAR) WW(II+1,JJ,KK) = WP
         IF (KK==0)    UU(II,JJ,KK)   = UM
         IF (KK==KBAR) UU(II,JJ,KK+1) = UP
         IF (CORRECTOR .AND. II>0 .AND. II<IBAR .AND. KK>0 .AND. KK<KBAR) THEN
            U_Z(II,JJ,KK) = 0.5_EB*(UM+UP) 
            W_X(II,JJ,KK) = 0.5_EB*(WM+WP) 
         ENDIF

      CASE(3) COMPONENT

         ICMM = CELL_INDEX(II,JJ,KK)
         ICMP = CELL_INDEX(II,JJ+1,KK)
         ICPM = CELL_INDEX(II+1,JJ,KK)
         ICPP = CELL_INDEX(II+1,JJ+1,KK)
         UP   = UU(II,JJ+1,KK)
         UM   = UU(II,JJ,KK)
         VP   = VV(II+1,JJ,KK)
         VM   = VV(II,JJ,KK)

         IF (NOM1==NM) THEN
            IWM = WALL_INDEX(ICMM, 2) 
            IWP = WALL_INDEX(ICPM, 2) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) UP = BC*UM
               IF (BC>1.5_EB) UP = FVT*SF%VEL_T(1)
            ENDIF
         ENDIF
         IF (NOM1==-NM) THEN
            IWM = WALL_INDEX(ICMP,-2) 
            IWP = WALL_INDEX(ICPP,-2) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) UM = BC*UP
               IF (BC>1.5_EB) UM = FVT*SF%VEL_T(1)
            ENDIF
         ENDIF

         IF (NOM2==NM) THEN
            IWM = WALL_INDEX(ICMM, 1) 
            IWP = WALL_INDEX(ICMP, 1) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) VP = BC*VM
               IF (BC>1.5_EB) VP = FVT*SF%VEL_T(1)
            ENDIF
         ENDIF
         IF (NOM2==-NM) THEN
            IWM = WALL_INDEX(ICPM,-1) 
            IWP = WALL_INDEX(ICPP,-1) 
            IF (BOUNDARY_TYPE(IWM)/=NULL_BOUNDARY .AND. BOUNDARY_TYPE(IWP)/=NULL_BOUNDARY) THEN
               IF (BC<1.5_EB) VM = BC*VP
               IF (BC>1.5_EB) VM = FVT*SF%VEL_T(1)
            ENDIF
         ENDIF

         IF (ABS(NOM1)/=NM .AND. NOM1/=0) THEN
            IF (NOM1<0) UM = OMESH(ABS(NOM1))%U(IIO1,JJO1,KKO1)
            IF (NOM1>0) UP = OMESH(ABS(NOM1))%U(IIO1,JJO1,KKO1)
         ENDIF
         IF (ABS(NOM2)/=NM .AND. NOM2/=0) THEN
            IF (NOM2<0) VM = OMESH(ABS(NOM2))%V(IIO2,JJO2,KKO2)
            IF (NOM2>0) VP = OMESH(ABS(NOM2))%V(IIO2,JJO2,KKO2)
         ENDIF

         MUA = .25_EB*( MU(II,JJ,KK) + MU(II+1,JJ,KK) + MU(II+1,JJ+1,KK) + MU(II,JJ+1,KK) )
         DVDX = RDXN(II)*(VP-VM)
         DUDY = RDYN(JJ)*(UP-UM)
         OME_E(IE) = DVDX - DUDY
         TAU_E(IE) = MUA*(DVDX + DUDY)
         IF (II==0)    VV(II,JJ,KK)   = VM
         IF (II==IBAR) VV(II+1,JJ,KK) = VP
         IF (JJ==0)    UU(II,JJ,KK)   = UM
         IF (JJ==JBAR) UU(II,JJ+1,KK) = UP
         IF (CORRECTOR .AND. II>0 .AND. II<IBAR .AND. JJ>0 .AND. JJ<JBAR) THEN
            U_Y(II,JJ,KK) = 0.5_EB*(UM+UP) 
            V_X(II,JJ,KK) = 0.5_EB*(VM+VP) 
         ENDIF

   END SELECT COMPONENT

   ! For SAWTOOTH=.FALSE., zero out vorticity and strain at boundary

   IF (IOR==0) THEN
      OME_E(IE) = 0._EB
      TAU_E(IE) = 0._EB
   ENDIF

ENDDO EDGE_LOOP

! Store cell node averages of the velocity components in UVW_GHOST for use in Smokeview

IF (CORRECTOR) THEN
   DO K=0,KBAR
      DO J=0,JBAR
         I_LOOP: DO I=0,IBAR
            IC = CELL_INDEX(I,J,K) 
            IF (IC==0) CYCLE I_LOOP
            IF (U_Y(I,J,K)  >-1.E5_EB) UVW_GHOST(IC,1) = U_Y(I,J,K) 
            IF (U_Y(I,J,K+1)>-1.E5_EB) UVW_GHOST(IC,1) = U_Y(I,J,K+1) 
            IF (U_Z(I,J,K)  >-1.E5_EB) UVW_GHOST(IC,1) = U_Z(I,J,K) 
            IF (U_Z(I,J+1,K)>-1.E5_EB) UVW_GHOST(IC,1) = U_Z(I,J+1,K) 
            IF (V_X(I,J,K)  >-1.E5_EB) UVW_GHOST(IC,2) = V_X(I,J,K) 
            IF (V_X(I,J,K+1)>-1.E5_EB) UVW_GHOST(IC,2) = V_X(I,J,K+1) 
            IF (V_Z(I,J,K)  >-1.E5_EB) UVW_GHOST(IC,2) = V_Z(I,J,K) 
            IF (V_Z(I+1,J,K)>-1.E5_EB) UVW_GHOST(IC,2) = V_Z(I+1,J,K) 
            IF (W_X(I,J,K)  >-1.E5_EB) UVW_GHOST(IC,3) = W_X(I,J,K) 
            IF (W_X(I,J+1,K)>-1.E5_EB) UVW_GHOST(IC,3) = W_X(I,J+1,K) 
            IF (W_Y(I,J,K)  >-1.E5_EB) UVW_GHOST(IC,3) = W_Y(I,J,K) 
            IF (W_Y(I+1,J,K)>-1.E5_EB) UVW_GHOST(IC,3) = W_Y(I+1,J,K) 
       !!!  UVW_GHOST(IC,1) = 0.25_EB*(U_Y(I,J,K)+U_Y(I,J,K+1)+U_Z(I,J,K)+U_Z(I,J+1,K)) 
       !!!  UVW_GHOST(IC,2) = 0.25_EB*(V_X(I,J,K)+V_X(I,J,K+1)+V_Z(I,J,K)+V_Z(I+1,J,K)) 
       !!!  UVW_GHOST(IC,3) = 0.25_EB*(W_X(I,J,K)+W_X(I,J+1,K)+W_Y(I,J,K)+W_Y(I+1,J,K)) 
         ENDDO I_LOOP
      ENDDO
   ENDDO
ENDIF

END SUBROUTINE VELOCITY_BC 
 
 
SUBROUTINE CHECK_STABILITY(NM)
 
! Checks the Courant and Von Neumann stability criteria, and if necessary, reduces the time step accordingly
 
REAL(EB) :: UODX,VODY,WODZ,UVW,UVWMAX,R_DX2,MU_MAX,MUTRM
INTEGER  :: NM,I,J,K
 
CHANGE_TIME_STEP(NM) = .FALSE.
UVWMAX = 0._EB
VN     = 0._EB
MUTRM  = 1.E-9_EB
R_DX2  = 1.E-9_EB
 
! Determine max CFL number from all grid cells
 
DO K=0,KBAR
   DO J=0,JBAR
      DO I=0,IBAR
         UODX = ABS(US(I,J,K))*RDXN(I)
         VODY = ABS(VS(I,J,K))*RDYN(J)
         WODZ = ABS(WS(I,J,K))*RDZN(K)
         UVW  = MAX(UODX,VODY,WODZ)
         IF (UVW>=UVWMAX) THEN
            UVWMAX = UVW 
            ICFL=I
            JCFL=J
            KCFL=K
         ENDIF
      ENDDO 
   ENDDO   
ENDDO   
 
CFL = DT*UVWMAX
 
! Determine max Von Neumann Number for fine grid calcs
 
PARABOLIC_IF: IF (DNS .OR. DXMIN<0.005_EB) THEN
 
   INCOMPRESSIBLE_IF: IF (ISOTHERMAL .AND. N_SPECIES==0) THEN
      IF (TWO_D) THEN
         R_DX2 = 1._EB/DXMIN**2 + 1._EB/DZMIN**2
      ELSE
         R_DX2 = 1._EB/DXMIN**2 + 1._EB/DYMIN**2 + 1._EB/DZMIN**2
      ENDIF
      MUTRM = MAX(RPR,RSC)*SPECIES(0)%MU(NINT(0.1_EB*TMPA))/RHOA
   ELSE INCOMPRESSIBLE_IF
      MU_MAX = 0._EB   
      DO K=1,KBAR
         DO J=1,JBAR
            IILOOP: DO I=1,IBAR
               IF (SOLID(CELL_INDEX(I,J,K))) CYCLE IILOOP
               IF (MU(I,J,K)>=MU_MAX) THEN
                  MU_MAX = MU(I,J,K)
                  I_VN=I
                  J_VN=J
                  K_VN=K
               ENDIF
            ENDDO IILOOP
         ENDDO  
      ENDDO  
      IF (TWO_D) THEN
         R_DX2 = RDX(I_VN)**2 + RDZ(K_VN)**2
      ELSE
         R_DX2 = RDX(I_VN)**2 + RDY(J_VN)**2 + RDZ(K_VN)**2
      ENDIF
      MUTRM = MAX(RPR,RSC)*MU_MAX/RHOS(I_VN,J_VN,K_VN)
   ENDIF INCOMPRESSIBLE_IF
 
   VN = DT*2._EB*R_DX2*MUTRM
 
ENDIF PARABOLIC_IF
 
! Adjust time step size if necessary
 
IF (CFL<CFL_MAX .AND. VN<VN_MAX) THEN
   DTNEXT = DT
   IF (CFL<=CFL_MIN .AND. VN<VN_MIN) DTNEXT = MIN(1.1_EB*DT,DTINT)
ELSE
   IF (UVWMAX==0._EB) UVWMAX = 1._EB
   DT = 0.9_EB*MIN( CFL_MAX/UVWMAX , VN_MAX/(2._EB*R_DX2*MUTRM) )
   CHANGE_TIME_STEP(NM) = .TRUE.
ENDIF
 
END SUBROUTINE CHECK_STABILITY
 
 
SUBROUTINE BAROCLINIC_CORRECTION
 
REAL(EB), POINTER, DIMENSION(:,:,:) :: UU,VV,WW,HQS,RTRM,RHOP
REAL(EB) :: RHO_AVG_OLD,RMIN,RMAX,RRAT,U2,V2,W2
INTEGER  :: I,J,K
 
HQS  => WORK1
RTRM => WORK2
 
IF (PREDICTOR) THEN
   UU => U
   VV => V
   WW => W
   RHOP=>RHO
ELSE
   UU => US
   VV => VS
   WW => WS
   RHOP=>RHOS
ENDIF
 
RHO_AVG_OLD = RHO_AVG
RMIN =  1000._EB
RMAX = -1000._EB
DO K=1,KBAR
      DO J=1,JBAR
         DO I=1,IBAR
            IF (.NOT.SOLID(CELL_INDEX(I,J,K))) THEN
               RMIN = MIN(RHOP(I,J,K),RMIN)
               RMAX = MAX(RHOP(I,J,K),RMAX)
            ENDIF
         ENDDO
      ENDDO
ENDDO
 
RHO_AVG = 2._EB*RMIN*RMAX/(RMIN+RMAX)
RRAT = RHO_AVG_OLD/RHO_AVG
RTRM = (1._EB-RHO_AVG/RHOP)*RRAT
 
DO K=1,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         U2 = 0.25_EB*(UU(I,J,K)+UU(I-1,J,K))**2
         V2 = 0.25_EB*(VV(I,J,K)+VV(I,J-1,K))**2
         W2 = 0.25_EB*(WW(I,J,K)+WW(I,J,K-1))**2
         HQS(I,J,K) = 0.5_EB*(U2+V2+W2)
      ENDDO
   ENDDO
ENDDO
 
DO K=1,KBAR
   DO J=1,JBAR
      U2 = (1.5_EB*UU(0,J,K)-0.5_EB*UU(1,J,K))**2
      V2 = 0.25_EB*(VV(1,J,K)+VV(1,J-1,K))**2
      W2 = 0.25_EB*(WW(1,J,K)+WW(1,J,K-1))**2
      HQS(0,J,K) = MIN(0.5_EB*(U2+V2+W2),10._EB*DX(1)+HQS(1,J,K))
      U2 = (1.5_EB*UU(IBAR,J,K)-0.5_EB*UU(IBM1,J,K))**2
      V2 = 0.25_EB*(VV(IBAR,J,K)+VV(IBAR,J-1,K))**2
      W2 = 0.25_EB*(WW(IBAR,J,K)+WW(IBAR,J,K-1))**2
      HQS(IBP1,J,K) = MIN(0.5_EB*(U2+V2+W2),10._EB*DX(IBAR)+HQS(IBAR,J,K))
   ENDDO
ENDDO
 
IF (.NOT.TWO_D) THEN
   DO K=1,KBAR
      DO I=1,IBAR
         U2 = 0.25_EB*(UU(I,1,K)+UU(I-1,1,K))**2
         V2 = (1.5_EB*VV(I,0,K)-0.5_EB*VV(I,1,K))**2
         W2 = 0.25_EB*(WW(I,1,K)+WW(I,1,K-1))**2
         HQS(I,0,K) = MIN(0.5_EB*(U2+V2+W2),10._EB*DY(1)+HQS(I,1,K))
         U2 = 0.25_EB*(UU(I,JBAR,K)+UU(I-1,JBAR,K))**2
         V2 = (1.5_EB*VV(I,JBAR,K)-0.5_EB*VV(I,JBM1,K))**2
         W2 = 0.25_EB*(WW(I,JBAR,K)+WW(I,JBAR,K-1))**2
         HQS(I,JBP1,K) = MIN(0.5_EB*(U2+V2+W2),10._EB*DY(JBAR)+HQS(I,JBAR,K))
      ENDDO
   ENDDO
ENDIF
 
DO J=1,JBAR
   DO I=1,IBAR
      U2 = 0.25_EB*(UU(I,J,1)+UU(I-1,J,1))**2
      V2 = 0.25_EB*(VV(I,J,1)+VV(I,J-1,1))**2
      W2 = (1.5_EB*WW(I,J,0)-0.5_EB*WW(I,J,1))**2
      HQS(I,J,0) = MIN(0.5_EB*(U2+V2+W2),10._EB*DZ(1)+HQS(I,J,1))
      U2 = 0.25_EB*(UU(I,J,KBAR)+UU(I-1,J,KBAR))**2
      V2 = 0.25_EB*(VV(I,J,KBAR)+VV(I,J-1,KBAR))**2
      W2 = (1.5_EB*WW(I,J,KBAR)-0.5_EB*WW(I,J,KBM1))**2
      HQS(I,J,KBP1) = MIN(0.5_EB*(U2+V2+W2),10._EB*DZ(KBAR)+HQS(I,J,KBAR))
   ENDDO
ENDDO
 
DO K=1,KBAR
   DO J=1,JBAR
      DO I=0,IBAR
         FVX(I,J,K) = FVX(I,J,K) - 0.5_EB*(RTRM(I+1,J,K)+RTRM(I,J,K))*(H(I+1,J,K)-H(I,J,K)-HQS(I+1,J,K)+HQS(I,J,K))*RDXN(I)
      ENDDO
   ENDDO
ENDDO
 
IF (.NOT.TWO_D) THEN
   DO K=1,KBAR
      DO J=0,JBAR
         DO I=1,IBAR
            FVY(I,J,K) = FVY(I,J,K) - 0.5_EB*(RTRM(I,J+1,K)+RTRM(I,J,K))*(H(I,J+1,K)-H(I,J,K)-HQS(I,J+1,K)+HQS(I,J,K))*RDYN(J)
         ENDDO
      ENDDO
   ENDDO
ENDIF
 
DO K=0,KBAR
   DO J=1,JBAR
      DO I=1,IBAR
         FVZ(I,J,K) = FVZ(I,J,K) - 0.5_EB*(RTRM(I,J,K+1)+RTRM(I,J,K))*(H(I,J,K+1)-H(I,J,K)-HQS(I,J,K+1)+HQS(I,J,K))*RDZN(K)
      ENDDO
   ENDDO
ENDDO
 
END SUBROUTINE BAROCLINIC_CORRECTION

SUBROUTINE GET_REV_velo(MODULE_REV,MODULE_DATE)
INTEGER,INTENT(INOUT) :: MODULE_REV
CHARACTER(255),INTENT(INOUT) :: MODULE_DATE

WRITE(MODULE_DATE,'(A)') velorev(INDEX(velorev,':')+1:LEN_TRIM(velorev)-2)
READ (MODULE_DATE,'(I5)') MODULE_REV
WRITE(MODULE_DATE,'(A)') velodate

END SUBROUTINE GET_REV_velo
 
END MODULE VELO