PROGRAM Vertical_Modes
 
USE m_Parameters

USE m_Arrays

USE m_Get_Vertical_Modes, ONLY: Create_Symmetric_Matrix, Get_Vertical_Modes

IMPLICIT NONE
 
! Computes the Vertical Mode Structure in Sigma Coordinate
! With the Equivalent Height as the Eigenvalue.
! The Vertical Discretization is Staggered.

! Input Variables:

! kmax  - Number of Vertical Levels
! ts    - Surface Temperature of Basic State (Kelvin)
! ps    - Surface Pressure of Basic State (mb)
! PrsLev(k) - Pressure Layers
! t0(k) - Basic State Temperature (Kelvin)
! gama(k) - Static Stability Parameter (m/hPa2

! Output Variables:

! dn(n)   - Equivalent Heights (Eigenvalues)
! fn(k,n) - Vertical Structure of the Modes (Eigenvectors)

! Others:

! z0(k)   - Hydrostatic Geopotential Height (From T0(K))
! z0s(1)  - Mean Surface Geopotentcial Height (From T0S(1))
! rd      - Dry Air Gas Constant
! gv      - Gravity Aceleration
! cp      - Specific Heat at Constant Pressure For Dry Air

INTEGER (KIND=ki8) :: l, n

REAL (KIND=kr4) :: bds

REAL (KIND=kr8) :: sinval, aux1, aux2, qs

REAL (KIND=kr8), DIMENSION (1) :: pp


CALL Get_Parameters

CALL Get_Arrays

OPEN (UNIT=10, FILE=TRIM(dirvmd)//dnfns, STATUS='UNKNOWN', FORM='UNFORMATTED')
bds=ps
WRITE (UNIT=10) bds
bds=ts
WRITE (UNIT=10) bds
bds=zs
WRITE (UNIT=10) bds
bdn=delprs
DO k=kmaxp,1,-1
  WRITE (UNIT=10) bdn(k)
END DO
bdn=prsint
DO k=kmaxp,1,-1
  WRITE (UNIT=10) bdn(k)
END DO
bdn=t0
DO k=kmaxp,1,-1
  WRITE (UNIT=10) bdn(k)
END DO
bdn=z0
DO k=kmaxp,1,-1
  WRITE (UNIT=10) bdn(k)
END DO
bdn=gama
DO k=kmaxp,1,-1
  WRITE (UNIT=10) bdn(k)
END DO
CLOSE (UNIT=10)
OPEN  (UNIT=15, FILE=TRIM(dirvmd)//dnfns//'.ctl', STATUS='UNKNOWN')
WRITE (UNIT=15, FMT='(A)') 'dset ^'//dnfns
WRITE (UNIT=15, FMT='(A)') 'undef 1.0e+20'
WRITE (UNIT=15, FMT='(A)') 'options sequential'
WRITE (UNIT=15, FMT='(A)') 'title Vertical Standard Structure'
WRITE (UNIT=15, FMT='(A)') 'xdef   1 linear 0 1'
WRITE (UNIT=15, FMT='(A)') 'ydef   1 linear 0 1'
WRITE (UNIT=15, FMT='(A,I3,A)') 'zdef ',kmaxp,' levels '
WRITE (UNIT=15, FMT='(7F9.3)') (prsint(k),k=kmaxp,1,-1)
WRITE (UNIT=15, FMT='(A)') 'tdef   1 linear 00z01jan2000 1dy'
WRITE (UNIT=15, FMT='(A)') 'vars 8'
WRITE (UNIT=15, FMT='(A)') 'ps    0 99 Standard Surface Pressure [hPa]'
WRITE (UNIT=15, FMT='(A)') 't0s   0 99 Standard Surface Temperature [K]'
WRITE (UNIT=15, FMT='(A)') 'z0s   0 99 Standard Surface Height [m]'
WRITE (UNIT=15, FMT='(A,I3,A)') 'del ',kmaxp,' 99 Delta Pressure [no dim]'
WRITE (UNIT=15, FMT='(A,I3,A)') 'prsint  ',kmaxp,' 99 Pressure of Layers [no dim]'
WRITE (UNIT=15, FMT='(A,I3,A)') 't0  ',kmaxp,' 99 Standard Temperature of Layers [K]'
WRITE (UNIT=15, FMT='(A,I3,A)') 'z0  ',kmaxp,' 99 Standard Height of Layers [m]'
WRITE (UNIT=15, FMT='(A,I3,A)') 'gama ',kmaxp,' 99 Static Stability Parameter [m/hPa2]'
WRITE (UNIT=15, FMT='(A)') 'endvars'
CLOSE (UNIT=15)

! Get the Matrix for Eigen-Problem

CALL Create_Symmetric_Matrix (mdiag,msubdiag)

WRITE (UNIT=*, FMT='(/,A)') ' mdiag:'
WRITE (UNIT=*, FMT='(1P6G12.5)') mdiag
WRITE (UNIT=*, FMT='(/,A)') ' msubdiag:'
WRITE (UNIT=*, FMT='(1P6G12.5)') msubdiag

! Eigenvalues and Eigenvectors

CALL Get_Vertical_Modes (mdiag, msubdiag, dn, fn)

! Results Ouput

OPEN (UNIT=20, FILE=TRIM(dirvmd)//dnfnb, STATUS='UNKNOWN', FORM='UNFORMATTED')
WRITE (UNIT=20) kmax
DO n=1,kmax
  WRITE (UNIT=20) dn(n)
  WRITE (UNIT=20) (fn(l,n),l=1,kmax)
END DO
WRITE (UNIT=20) PrsLev(1:kmax)
WRITE (UNIT=20) z0(1:kmax)
CLOSE (UNIT=20)

OPEN (UNIT=30, FILE=TRIM(dirvmd)//dnfng, STATUS='UNKNOWN', FORM='UNFORMATTED')
DO n=1,kmax
  bdm(n)=dn(n)
END DO
WRITE (UNIT=30) bdm
DO k=kmax,1,-1
  DO n=1,kmax
    bdm(n)=fn(k,n)
  END DO
  WRITE (UNIT=30) bdm
END DO
CLOSE (UNIT=30)

OPEN  (UNIT=35, FILE=TRIM(dirvmd)//dnfng//'.ctl', STATUS='UNKNOWN')
WRITE (UNIT=35, FMT='(A)') 'dset ^'//dnfng
WRITE (UNIT=35, FMT='(A)') 'undef 1.0e+20'
WRITE (UNIT=35, FMT='(A)') 'options sequential'
WRITE (UNIT=35, FMT='(A)') 'title Vertical Modes Structure'
WRITE (UNIT=35, FMT='(A,I3,A)') 'xdef ',kmax,' linear 0 1'
WRITE (UNIT=35, FMT='(A)') 'ydef   1 linear 0 1'
WRITE (UNIT=35, FMT='(A,I3,A)') 'zdef ',kmax,' levels '
WRITE (UNIT=35, FMT='(7F9.3)') (PrsLev(k),k=kmax,1,-1)
WRITE (UNIT=35, FMT='(A)') 'tdef   1 linear 00z01jan2000 1dy'
WRITE (UNIT=35, FMT='(A)') 'vars 2'
WRITE (UNIT=35, FMT='(A)') 'dn   0 99 Equivalent Depth [m]'
WRITE (UNIT=35, FMT='(A,I3,A)') 'fn ',kmax,' 99 Vertical Eigenfunction'
WRITE (UNIT=35, FMT='(A)') 'endvars'
CLOSE (UNIT=35)

WRITE (UNIT=*, FMT='(A)') ' '
WRITE (UNIT=*, FMT='(A)') ' Ortogonalidade'
WRITE (UNIT=*, FMT='(A)') ' '
DO n=1,kmax
  DO l=1,kmax
    or(n,l)=0.0
    DO k=1,kmax
      or(n,l)=or(n,l)+fn(k,n)*fn(k,l)
    END DO
    IF (n == l) THEN
      IF (ABS(1.0_kr8-or(n,l)) < 1.0E-12_kr8) or(n,l)=1.0_kr8
    ELSE 
      IF (ABS(or(n,l)) < 1.0E-12_kr8) or(n,l)=0.0_kr8
    END IF
    IF (or(n,l) /= 0.0_kr8) WRITE (UNIT=*, FMT='(2I5,1PG16.8)') n, l, or(n,l)
  END DO
END DO

DEALLOCATE (PrsLev, prsint, delprs, t0, z0)
CALL Free_Arrays


END PROGRAM Vertical_Modes