PROGRAM BLGwind

        IMPLICIT NONE

        INTEGER n,i,j,k,imax,jmax,kmax
        REAL z, Cd, fc, delta_T, Ui, Vi, delta_T1,VMAX,CDIN,UIMB,VIMB,UDIFF,VDIFF
        REAL, Allocatable :: U(:,:,:), V(:,:,:), M(:,:,:), delta_U(:,:,:), delta_V(:,:,:),Vgr(:,:,:),R(:)

        OPEN(6,file='output.dat',status='unknown',position='append')
        OPEN(5,file='input.dat',status='old')
        READ(5,*)imax,jmax,kmax,VMAX,Z
        WRITE(6,*)
        WRITE(6,37)VMAX,Z
        WRITE(6,*)

        Allocate(R(jmax))
        Allocate(U(imax,jmax,kmax), V(imax,jmax,kmax), M(imax,jmax,kmax))
        Allocate(delta_U(imax,jmax,kmax), delta_V(imax,jmax,kmax),Vgr(imax,jmax,kmax))

!       z       = 2000. !1500.0     !Height of Atmosphere(m)
!       WRITE(0,*)'INPUT Zi in m'
!       READ(0,*)z
!       Cd      = CDIN !0.001 + (4.0/100000.)*VMAX !0.0020 !CDIN !0.0020       !(drag coefficient)
!
        fc      = 0.00015    !Coriolis Parameter
        delta_T1= 60.0     !time step(sec)(1 hr)

        DO j=1,jmax         ! j along the radius, i along azimuth
          R(j) = 90000. !60000.0   !Radius of Curvature(m)
        ENDDO

        DO k=1,kmax
         DO j=1,jmax
          DO i=1,imax
            Vgr(i,j,k)= VMAX                ! tangential wind speed at the top of the PBL 
          ENDDO
         ENDDO
        ENDDO

        WRITE(*,*) " "
        WRITE(*,*) "   BOUNDARY LAYER GRADIENT WIND!"
        WRITE(*,*) "Iteration      Vt(m/s)     Ur(m/s)     M(m/s)   delta U     delta V"

!       Assume U and V equal zero first

        n=0
        Ui=0
        Vi=VMAX

        DO k=1,kmax
         DO j=1,jmax
          DO i=1,imax
           M(i,j,k)=Vi
          ENDDO
         ENDDO
        ENDDO
        DO k=1,kmax
         DO j=1,jmax
          DO i=1,imax 
             DO WHILE (n.LE.300)     ! iterative loop

                IF(M(i,j,k) .LE. 20.0)THEN
                  Cd= 2.0*(0.7/1000.) + (6.5/100000.)*M(i,j,k) !VMAX
                ELSE
                  Cd=2.0*0.0020
                ENDIF

                IF (n.EQ.0) THEN
                        delta_T = 0.0
                ELSE IF (n>0) THEN
                        delta_T =delta_T1 
                ENDIF

!               UIMB     = (-Vgr(i,j,k) + (fc*Vi) + (Vi*Vi/R(j)))
                UIMB     = -(Vgr(i,j,k)**2 - Vi**2)/R(j) - fc*(Vgr(i,j,k) - Vi)
                UDIFF    = (Cd*M(i,j,k)*Ui/z)
                U(i,j,k) = Ui + delta_T * (UIMB - UDIFF)

!               VIMB     = (-(fc*Ui)            - (Ui*Vi/R(j)))
                VIMB     = -(Vi/R(j) + fc)*Ui 
                VDIFF    = (Cd*M(i,j,k)*Vi/z)
                V(i,j,k) = Vi + delta_T * (VIMB - VDIFF)

                IF(n.EQ.0)THEN  
                        delta_U = 0.0
                ELSE IF(n>0)THEN
                        delta_U(i,j,k) = U(i,j,k) - Ui
                        delta_V(i,j,k) = V(i,j,k) - Vi
                ENDIF

                M(i,j,k) = SQRT((U(i,j,k))**2+(V(i,j,k))**2)

                Ui = U(i,j,k)
                Vi = V(i,j,k)


                n=n+1

                if(i==1 .and. j==1 .and. k==1)WRITE(*,38) n,  V(i,j,k), -U(i,j,k), M(i,j,k), delta_U(i,j,k), delta_V(i,j,k)
                if(i==1 .and. j==1 .and. k==1)WRITE(6,39)n,-U(i,j,k),V(i,j,k) !,-3600.0*UIMB,3600.0*UDIFF,3600.0*VIMB,3600.0*VDIFF
             ENDDO
          ENDDO
         ENDDO
        ENDDO

        WRITE(6,*)
37      FORMAT(1X,F4.1,4X,F8.1)
38      FORMAT(4X,I4,10X,F6.2,6X,F6.2,6X,F5.2,6X,F5.3,5X,F6.3)
39      FORMAT(I5,4X,6(F8.4,2x))
        DeAllocate (U, V, M, delta_U, delta_V,Vgr,R)

        END