PROGRAM Vertical_Decomposition

USE m_Parameters

USE m_Arrays

IMPLICIT NONE

! Perfoms Vertical Decomposition of
! Wind and Geopotential Height in Sigma Level.

! Input Variables:

! Imax    - Number of Longitudes
! Jmax    - Number of Latitudes
! Kmax    - Number of Vertical Levels

! uobs(i,j,k) - Observed Zonal Wind (m/s)
! vobs(i,j,k) - Observed Meridional Wind (m/s)
! hobs(i,j,k) - Observed Geopotential Height Deviation (m) [zobs - z0]
! z0(k)       - Basic State Geopotential Height (m)
! dn(n)       - Equivalent Heights (Eigenvalues)
! fn(k,n)     - Vertical Structure of the Modes (Eigenvectors)

! Output Variables:

! un(i,j,n) - Projection of Zonal Wind
! vn(i,j,n) - Projection of Meridional Wind
! hn(i,j,n) - Projection of Modified Geopotential Height
! ec(i,j,n) - Kinetic Energy of Each Mode
! ep(i,j,n) - Available Potential Energy of Each Mode
! et(i,j,n) - Total Energy of Each Mode
! pe(i,j,n) - Percentage of Energy of Each Mode
! ect(i,j)  - Total Kinetic Energy
! ept(i,j)  - Total Available Potential Energy
! ett(i,j)  - Total Energy

! Others:

! rd - Dry Air Gas Constant
! gv - Gravity Aceleration

INTEGER, PARAMETER :: nvarp=15, nvarn=23

INTEGER (KIND=ki4):: nr, lrec, &
                     nc0a, nc0b, nc1a, nc1b, nc2a, nc2b, &
                     nc3a, nc3b, nc4a, nc4b

INTEGER (KIND=ki8):: k, l, i, j, kk, n, m, Nmax

CHARACTER (LEN=5) :: cls

REAL (KIND=kr4) :: bgm
REAL (KIND=kr8) :: rad, p2, aux1, aux2, dx, etnn, etnm

CHARACTER (LEN=4) :: etcn(nvarn)
CHARACTER (LEN=5) :: etcp(nvarp)
CHARACTER (LEN=7) :: ncls(nvarn)

CALL Get_Parameters

CALL Get_Arrays

! Vertical Modes Input

OPEN (UNIT=10, FILE=TRIM(dirvmd)//dnfnb, STATUS='UNKNOWN', FORM='UNFORMATTED')
READ (UNIT=10) Nmax
IF (Nmax /= Kmax) THEN
  WRITE (UNIT=*, FMT='(A,I5)') ' The Number of Vertical Modes  = ',Nmax
  WRITE (UNIT=*, FMT='(A)')    ' Is Different than '
  WRITE (UNIT=*, FMT='(A,I5)') ' The Number of Vertical Layers = ',Kmax
  STOP 99
END IF
DO n=1,Kmax
  READ (UNIT=10) dn(n)
  cn2(n)=gv*dn(n)
  cn(n)=SQRT(cn2(n))
  READ (UNIT=10) (fn(k,n),k=1,Kmax)
END DO
READ  (UNIT=10) prslev
READ  (UNIT=10) z0
CLOSE (UNIT=10)

CALL Get_Classes

! Data Input

OPEN (UNIT=20, FILE=TRIM(dirdpr)//datain, STATUS='OLD', FORM= 'UNFORMATTED')
DO k=Kmax,1,-1
  READ (UNIT=20) buf
  DO j=1,Jmax
    DO i=1,Imax
      uobs(i,j,k)=REAL(buf(i,j),kr8)
    END DO
  END DO
END DO
DO k=Kmax,1,-1
  READ (UNIT=20) buf
  DO j=1,Jmax
    DO i=1,Imax
      vobs(i,j,k)=REAL(buf(i,j),kr8)
    END DO
  END DO
END DO
DO k=Kmax,1,-1
  READ (UNIT=20) buf
  DO j=1,Jmax
    DO i=1,Imax
      hobs(i,j,k)=REAL(buf(i,j),kr8)-z0(k)
    END DO
  END DO
END DO
CLOSE (UNIT=20)

! Fourier Latitudes and Gaussian Latitudes

dx=360.0/REAL(Imax,kr8)
CALL Get_Gaussian_Latitudes

! Data Do-Loop

DO i=1,Imax
  DO j=1,Jmax

    ! Vertical Decomposition

    DO n=1,Kmax
      un(i,j,n)=0.0
      vn(i,j,n)=0.0
      hn(i,j,n)=0.0
      DO k=1,Kmax
        aux1=fn(k,n)/cn(n)
        un(i,j,n)=un(i,j,n)+uobs(i,j,k)*aux1
        vn(i,j,n)=vn(i,j,n)+vobs(i,j,k)*aux1
        hn(i,j,n)=hn(i,j,n)+hobs(i,j,k)*fn(k,n)/dn(n)
      END DO
    END DO

    ! Kinetic, Available Potential and Total Energy

    DO n=1,Kmax
      aux2=0.5*cn2(n)
      ec(i,j,n)=aux2*(un(i,j,n)*un(i,j,n)+vn(i,j,n)*vn(i,j,n))
      ep(i,j,n)=aux2*hn(i,j,n)*hn(i,j,n)
      et(i,j,n)=ec(i,j,n)+ep(i,j,n)
    END DO
    ect(i,j)=0.0
    ept(i,j)=0.0
    DO n=1,Kmax
      ect(i,j)=ect(i,j)+ec(i,j,n)
      ept(i,j)=ept(i,j)+ep(i,j,n)
    END DO
    ett(i,j)=ect(i,j)+ept(i,j)

    ! Percentage of Total Energy

    aux2=100.0/ett(i,j)
    DO n=1,Kmax
      pe(i,j,n)=et(i,j,n)*aux2
    END DO

  END DO
END DO

INQUIRE (IOLENGTH=lrec) buf(1:Imax,1:Jmax)
WRITE (UNIT=*, FMT='(/,A,I7,/)') 'lrec = ', lrec
WRITE (UNIT=0, FMT='(/,A,I7,/)') 'lrec = ', lrec
OPEN (UNIT=30, FILE=TRIM(dirvdo)//vrtdec, STATUS='UNKNOWN', &
      FORM='UNFORMATTED', ACCESS='DIRECT', RECL=lrec)
nr=0
DO n=1,Kmax
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(un(i,j,n),kr4)
    END DO
  END DO
  nr=nr+1
  WRITE (UNIT=30, REC=nr) buf
END DO
DO n=1,Kmax
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(vn(i,j,n),kr4)
    END DO
  END DO
  nr=nr+1
  WRITE (UNIT=30, REC=nr) buf
END DO
DO n=1,Kmax
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(hn(i,j,n),kr4)
    END DO
  END DO
  nr=nr+1
  WRITE (UNIT=30, REC=nr) buf
END DO
CLOSE (UNIT=30)

OPEN  (UNIT=40, FILE=TRIM(dirvdo)//vrtdecb, STATUS='UNKNOWN')
WRITE (UNIT=40, FMT='(A)') 'dset ^'//vrtdec
WRITE (UNIT=40, FMT='(A)') 'undef 1.0E+20'
WRITE (UNIT=40, FMT='(A,I3,A)') 'title Vertical Modes Decomposition  Lev = Hn'
WRITE (UNIT=40, FMT='(A,I4,A,F9.6)') 'xdef ', Imax, ' linear 0.0 ', dx
WRITE (UNIT=40, FMT='(A,I4,A)') 'ydef ', Jmax, ' levels '
WRITE (UNIT=40, FMT='(10F8.3)') (glat(j),j=1,Jmax)
WRITE (UNIT=40, FMT='(A,I3,A)') 'zdef ', Kmax, ' levels '
WRITE (UNIT=40, FMT='(6F12.5)') (dn(n),n=1,Kmax)
WRITE (UNIT=40, FMT='(A)') 'tdef 1 linear '//grdtim//' 1dy'
WRITE (UNIT=40, FMT='(A)') 'vars 3'
WRITE (UNIT=40, FMT='(A,I3,A)') 'un ', Kmax, ' 99 Zonal Wind [no Dim]'
WRITE (UNIT=40, FMT='(A,I3,A)') 'vn ', Kmax, ' 99 Meridional Wind [no Dim]'
WRITE (UNIT=40, FMT='(A,I3,A)') 'hn ', Kmax, ' 99 Modified Geopotential Height [no Dim]'
WRITE (UNIT=40, FMT='(A)') 'endvars'
CLOSE (UNIT=40)

OPEN (UNIT=50, FILE=TRIM(dirvdo)//energn, STATUS='UNKNOWN', FORM='UNFORMATTED')

buf(:,:)=REAL(ect(:,:),kr4)
WRITE (UNIT=50) buf
etcn(1)='ect '
ncls(1)='       '
buf(:,:)=REAL(ept(:,:),kr4)
WRITE (UNIT=50) buf
etcn(2)='ept '
ncls(2)='       '
buf(:,:)=REAL(ett(:,:),kr4)
WRITE (UNIT=50) buf
etcn(3)='ett '
ncls(3)='       '

CALL CnEc (nc0a, nc0b)
etcn(4)='ecn0'
ncls(4)='Class 0'
CALL CnEc (nc1a, nc1b)
etcn(5)='ecn1'
ncls(5)='Class 1'
CALL CnEc (nc2a, nc2b)
etcn(6)='ecn2'
ncls(6)='Class 2'
CALL CnEc (nc3a, nc3b)
etcn(7)='ecn3'
ncls(7)='Class 3'
CALL CnEc (nc4a, nc4b)
etcn(8)='ecn4'
ncls(8)='Class 4'

CALL CnEp (nc0a, nc0b)
etcn(9)='epn0'
ncls(9)='Class 0'
CALL CnEp (nc1a, nc1b)
etcn(10)='epn1'
ncls(10)='Class 1'
CALL CnEp (nc2a, nc2b)
etcn(11)='epn2'
ncls(11)='Class 2'
CALL CnEp (nc3a, nc3b)
etcn(12)='epn3'
ncls(12)='Class 3'
CALL CnEp (nc4a, nc4b)
etcn(13)='epn4'
ncls(13)='Class 4'

CALL CnEt (nc0a, nc0b)
etcn(14)='etn0'
ncls(14)='Class 0'
CALL CnEt (nc1a, nc1b)
etcn(15)='etn1'
ncls(15)='Class 1'
CALL CnEt (nc2a, nc2b)
etcn(16)='etn2'
ncls(16)='Class 2'
CALL CnEt (nc3a, nc3b)
etcn(17)='etn3'
ncls(17)='Class 3'
CALL CnEt (nc4a, nc4b)
etcn(18)='etn4'
ncls(18)='Class 4'

CALL CnPe (nc0a, nc0b)
etcn(19)='pen0'
ncls(19)='Class 0'
CALL CnPe (nc1a, nc1b)
etcn(20)='pen1'
ncls(20)='Class 1'
CALL CnPe (nc2a, nc2b)
etcn(21)='pen2'
ncls(21)='Class 2'
CALL CnPe (nc3a, nc3b)
etcn(22)='pen3'
ncls(22)='Class 3'
CALL CnPe (nc4a, nc4b)
etcn(23)='pen4'
ncls(23)='Class 4'

CLOSE (UNIT=50)

OPEN  (UNIT=60, FILE=TRIM(dirvdo)//energnb, STATUS='UNKNOWN')
WRITE (UNIT=60, FMT='(A)') 'dset ^'//energn
WRITE (UNIT=60, FMT='(A)') 'options sequential'
WRITE (UNIT=60, FMT='(A)') 'undef 1.0E+20'
WRITE (UNIT=60, FMT='(A)') 'title Vertical Modes Decomposition'
WRITE (UNIT=60, FMT='(A,I4,A,F9.6)') 'xdef ', Imax, ' linear 0.0 ', dx
WRITE (UNIT=60, FMT='(A,I4,A)') 'ydef ', Jmax, ' levels '
WRITE (UNIT=60, FMT='(10F8.3)') (glat(j),j=1,Jmax)
WRITE (UNIT=60, FMT='(A)') 'zdef 1 levels 1000'
WRITE (UNIT=60, FMT='(A)') 'tdef 1 linear '//grdtim//' 1dy'
WRITE (UNIT=60, FMT='(A,I3)') 'vars ', nvarn
WRITE (UNIT=60, FMT='(A)') etcn(1)//' 0 99 Total Kinetic Energy [m2/s2]'
WRITE (UNIT=60, FMT='(A)') etcn(2)//' 0 99 Total Potential Energy [m2/s2]'
WRITE (UNIT=60, FMT='(A)') etcn(3)//' 0 99 Total Total Energy [m2/s2]'
DO n=4,8
WRITE (UNIT=60, FMT='(A)') etcn(n)//' 0 99 Modal Kinetic Energy '//ncls(n)//' [m2/s2]'
END DO
DO n=9,13
WRITE (UNIT=60, FMT='(A)') etcn(n)//' 0 99 Modal Potential Energy '//ncls(n)//' [m2/s2]'
END DO
DO n=14,18
WRITE (UNIT=60, FMT='(A)') etcn(n)//' 0 99 Modal Total Energy '//ncls(n)// '[m2/s2]'
END DO
DO n=19,nvarn
WRITE (UNIT=60, FMT='(A)') etcn(n)//' 0 99 Modal Percentage of Energy '//ncls(n)//' [%]'
END DO
WRITE (UNIT=60, FMT='(A)') 'endvars'
CLOSE (UNIT=60)

OPEN (UNIT=70,FILE=TRIM(dirvdo)//enercl, STATUS='UNKNOWN', FORM='UNFORMATTED')
OPEN (UNIT=75,FILE=TRIM(dirvdo)//enercg, STATUS='UNKNOWN', FORM='UNFORMATTED')

!Class 0
CALL Cm (nc0a, nc0b)
etcp(1)='etc0 '

!Class 1
CALL Cm (nc1a, nc1b)
etcp(2)='etc1 '

!Class 2
CALL Cm (nc2a, nc2b)
etcp(3)='etc2 '

!Class 3
CALL Cm (nc3a, nc3b)
etcp(4)='etc3 '

!Class 4
CALL Cm (nc4a, nc4b)
etcp(5)='etc4 '

!classe 0 com 1
CALL CnCm (nc0a, nc0b, nc1a, nc1b)
etcp(6)='etc01'

!classe 0 com 2
CALL CnCm (nc0a, nc0b, nc2a, nc2b)
etcp(7)='etc02'

!classe 0 com 3
CALL CnCm (nc0a, nc0b, nc3a, nc3b)
etcp(8)='etc03'

!classe 0 com 4
CALL CnCm (nc0a, nc0b, nc4a, nc4b)
etcp(9)='etc04'

!classe 1 com 2
CALL CnCm (nc1a, nc1b, nc2a, nc2b)
etcp(10)='etc12'

!classe 1 com 3
CALL CnCm (nc1a, nc1b, nc3a, nc3b)
etcp(11)='etc13'

!classe 1 com 4
CALL CnCm (nc1a, nc1b, nc4a, nc4b)
etcp(12)='etc14'

!classe 2 com 3
CALL CnCm (nc2a, nc2b, nc3a, nc3b)
etcp(13)='etc23'

!classe 2 com 4
CALL CnCm (nc2a, nc2b, nc4a, nc4b)
etcp(14)='etc24'

!classe 3 com 4
CALL CnCm (nc3a, nc3b, nc4a, nc4b)
etcp(15)='etc34'

CLOSE (UNIT=70)
CLOSE (UNIT=75)

WRITE (UNIT=*, FMT='(/,A,I5,/)') ' Number of Variables for Interaction = ', nvarp
WRITE (UNIT=0, FMT='(/,A,I5,/)') ' Number of Variables for Interaction = ', nvarp

OPEN  (UNIT=80, FILE=TRIM(dirvdo)//enercgb, STATUS='UNKNOWN')
WRITE (UNIT=80, FMT='(A)') 'dset ^'//enercg
WRITE (UNIT=80, FMT='(A)') 'options sequential'
WRITE (UNIT=80, FMT='(A)') 'undef 1.0E+20'
WRITE (UNIT=80, FMT='(A)') 'title Vertical Modes Decomposition  Lev = Prs'
WRITE (UNIT=80, FMT='(A)') 'xdef 1 levels 0'
WRITE (UNIT=80, FMT='(A)') 'ydef 1 levels 0'
WRITE (UNIT=80, FMT='(A,I3,A)') 'zdef ', Kmax, ' levels '
WRITE (UNIT=80, FMT='(6F12.5)') (prslev(k),k=Kmax,1,-1)
WRITE (UNIT=80, FMT='(A)') 'tdef 1 linear '//grdtim//' 1dy'
WRITE (UNIT=80, FMT='(A,I5)') 'vars ', nvarp
DO n=1,5
  WRITE(UNIT=80, FMT='(2A,I3,A,I2,A)') etcp(n), ' ', Kmax,  &
                  ' 99 Modal Total Energy by Class', n-1, ' [m2/s2]'
END DO
DO n=6,nvarp
  cls=etcp(n)(4:4)//' - '//etcp(n)(5:5)
  WRITE (UNIT=80, FMT='(2A,I3,3A)') etcp(n), ' ', Kmax,  &
                     ' 99 Modal Total Energy - Interaction by Classes ', cls, ' [m2/s2]'
END DO
WRITE (UNIT=80, FMT='(A)') 'endvars'
CLOSE (UNIT=80)

OPEN  (UNIT=85, FILE=TRIM(dirvdo)//enerclb, STATUS='UNKNOWN')
WRITE (UNIT=85, FMT='(A)') 'dset ^'//enercl
WRITE (UNIT=85, FMT='(A)') 'options sequential'
WRITE (UNIT=85, FMT='(A)') 'undef 1.0E+20'
WRITE (UNIT=85, FMT='(A,I3,A)') 'title Vertical Modes Decomposition  Lev = Prs'
WRITE (UNIT=85, FMT='(A,I4,A,F9.6)') 'xdef ', Imax, ' linear 0.0 ', dx
WRITE (UNIT=85, FMT='(A,I4,A)') 'ydef ', Jmax, ' levels '
WRITE (UNIT=85, FMT='(10F8.3)') (glat(j),j=1,Jmax)
WRITE (UNIT=85, FMT='(A,I3,A)') 'zdef ', Kmax, ' levels '
WRITE (UNIT=85, FMT='(6F12.5)') (prslev(k),k=Kmax,1,-1)
WRITE (UNIT=85, FMT='(A)') 'tdef 1 linear '//grdtim//' 1dy'
WRITE (UNIT=85, FMT='(A,I5)') 'vars ', nvarp
DO n=1,5
  WRITE(UNIT=85, FMT='(2A,I3,A,I2,A)') etcp(n), ' ', Kmax,  &
                  ' 99 Modal Total Energy by Class', n-1, ' [m2/s2]'
END DO
DO n=6,nvarp
  cls=etcp(n)(4:4)//' - '//etcp(n)(5:5)
  WRITE (UNIT=85, FMT='(2A,I3,3A)') etcp(n), ' ', Kmax,  &
                     ' 99 Modal Total Energy - Interaction by Classes ', cls, ' [m2/s2]'
END DO
WRITE (UNIT=85, FMT='(A)') 'endvars'
CLOSE (UNIT=85)

CALL Free_Arrays


CONTAINS


SUBROUTINE Get_Gaussian_Latitudes

! Calculates Gaussian Latitudes

  IMPLICIT NONE

  INTEGER (KIND=ki8):: j

  REAL (KIND=kr8) :: eps, rd, dradz, drad, wgt, p1
  
  ! Gaussian Latitudes

  if (KIND(eps) .eq. kr8) then
     eps=1.0e-12_kr8
  else
     eps=1.0e-06
  endif
  rd=45.0/ATAN(1.0)
  dradz=((180.0/REAL(Jmax,kr8))/rd)/10.0
  wgt=2.0/(REAL(Jmax,kr8)*REAL(Jmax,kr8))
  rad=0.0
  DO j=1,Jmax/2
    drad=dradz
    DO WHILE (drad > eps)
      CALL Legendre_Polynomial (Jmax, rad, p2)
      DO
        p1=p2
        rad=rad+drad
        CALL Legendre_Polynomial (Jmax, rad, p2)
        IF (SIGN(1.0_kr8,p1) /= SIGN(1.0_kr8,p2)) EXIT
      END DO
      rad=rad-drad
      drad=drad*0.25
    END DO
    glat(j)=rad*rd-90.0
    glat(Jmax-j+1)=-glat(j)
    CALL Legendre_Polynomial (Jmax-1, rad, p1)
    wglat(j)=wgt*(1.0-COS(rad)*COS(rad))/(p1*p1)
    wglat(Jmax-j+1)=wglat(j)
    rad=rad+dradz
  END DO

END SUBROUTINE Get_Gaussian_Latitudes


SUBROUTINE Legendre_Polynomial (jdim, rad, p)

  ! Calculates the Value of the Ordinary Legendre Function
  ! at a Given Order Jmax and at a Specified Co-Latitude rad
  ! The Result Will Return in p

  IMPLICIT NONE

  INTEGER (KIND=ki8), INTENT (IN) :: jdim
  REAL (KIND=kr8), INTENT(IN) :: rad
  REAL (KIND=kr8), INTENT(OUT) :: p

  INTEGER (KIND=ki8):: n

  REAL (KIND=kr8) :: x, y1, y2, g, y3

  x=COS(rad)
  y1=1.0
  y2=X
  DO n=2,jdim
    g=x*y2
    y3=g-y1+g-(g-y1)/REAL(n,kr8)
    y1=y2
    y2=y3
  END DO
  p=y3

END SUBROUTINE Legendre_Polynomial


SUBROUTINE Global_Mean (buf,bgm)

  IMPLICIT NONE

  REAL (KIND=kr4), INTENT(OUT) :: bgm
  REAL (KIND=kr4), INTENT(IN) :: buf(Imax,Jmax)

  INTEGER :: j
  REAL (KIND=kr4) :: bzm, swgt

  swgt=SUM(wglat)
  bgm=0.0
  DO j=1,Jmax
    bzm=SUM(buf(:,j))/REAL(Imax,kr8)
    bgm=bgm+bzm*wglat(j)/swgt
  ENDDO

END SUBROUTINE Global_Mean


SUBROUTINE Get_Classes

  IMPLICIT NONE

  INTEGER (KIND=ki4) :: nc0, nc1, nc2, nc3, nc4, n
  REAL (KIND=kr8) :: Hc0, hC1a, Hc1b, Hc2a, Hc2b, Hc3a, Hc3b, Hc4

  ! Classes:
  ! Class 0 : Hn >= 600 m (Hc0)
  ! Class 1 : 100 m (Hc1a) =< Hn < 600 m (Hc1b=Hc0)
  ! Class 2 :  10 m (Hc2a) =< Hn < 100 m (Hc2b=Hc1a)
  ! Class 3 :   1 m (Hc3a) =< Hn <  10 m (Hc3b=Hc2a)
  ! Class 4 : Hn < 1 m (Hc4=Hc3a)

  Hc0=600.0
  Hc1a=100.0
  Hc1b=Hc0
  Hc2a=10.0
  Hc2b=Hc1a
  Hc3a=1.0
  Hc3b=Hc2a
  Hc4=Hc3a

  nc0=0
  nc1=0
  nc2=0
  nc3=0
  nc4=0

  DO n=1,Kmax
    ! Class 0
    IF (dn(n) >= Hc0) THEN
       nc0=nc0+1
    ENDIF
    ! Class 1
    IF (dn(n) >= Hc1a .AND. dn(n) < Hc1b) THEN
       nc1=nc1+1
    ENDIF
    ! Class 2
    IF (dn(n) >= Hc2a .AND. dn(n) < Hc2b) THEN
       nc2=nc2+1
    ENDIF
    ! Class 3
    IF (dn(n) >= Hc3a .AND. dn(n) < Hc3b) THEN
       nc3=nc3+1
    ENDIF
    ! Class 4
    IF (dn(n) < Hc4) THEN
       nc4=nc4+1
    ENDIF
  END DO

  nc0a=1
  nc0b=nc0
  nc1a=nc0b+1
  nc1b=nc1a+nc1-1
  nc2a=nc1b+1
  nc2b=nc2a+nc2-1
  nc3a=nc2b+1
  nc3b=nc3a+nc3-1
  nc4a=nc3b+1
  nc4b=nc4a+nc4-1

  WRITE (UNIT=*, FMT='(A)') ' '
  WRITE (UNIT=*, FMT='(3A,I5)') ' Caso = ', caso, ' Kmax = ', Kmax
  WRITE (UNIT=*, FMT='(A)') ' '
  WRITE (UNIT=*, FMT='(2(A,I5),A)') ' nc0a = ', nc0a, ' nc0b = ', nc0b, '  Classe 0'
  WRITE (UNIT=*, FMT='(2(A,I2,A,F9.4))') ' Hn(', nc0a,') = ', dn(nc0a), '  Hn(', nc0b,') = ', dn(nc0b)
  WRITE (UNIT=*, FMT='(A)') ' '
  WRITE (UNIT=*, FMT='(2(A,I5),A)') ' nc1a = ', nc1a, ' nc1b = ', nc1b, '  Classe 1'
  WRITE (UNIT=*, FMT='(2(A,I2,A,F9.4))') ' Hn(', nc1a,') = ', dn(nc1a), '  Hn(', nc1b,') = ', dn(nc1b)
  WRITE (UNIT=*, FMT='(A)') ' '
  WRITE (UNIT=*, FMT='(2(A,I5),A)') ' nc2a = ', nc2a, ' nc2b = ', nc2b, '  Classe 2'
  WRITE (UNIT=*, FMT='(2(A,I2,A,F9.4))') ' Hn(', nc2a,') = ', dn(nc2a), '  Hn(', nc2b,') = ', dn(nc2b)
  WRITE (UNIT=*, FMT='(A)') ' '
  WRITE (UNIT=*, FMT='(2(A,I5),A)') ' nc3a = ', nc3a, ' nc3b = ', nc3b, '  Classe 3'
  WRITE (UNIT=*, FMT='(2(A,I2,A,F9.4))') ' Hn(', nc3a,') = ', dn(nc3a), '  Hn(', nc3b,') = ', dn(nc3b)
  WRITE (UNIT=*, FMT='(A)') ' '
  WRITE (UNIT=*, FMT='(2(A,I5),A)') ' nc4a = ', nc4a, ' nc4b = ', nc4b, '  Classe 4'
  WRITE (UNIT=*, FMT='(2(A,I2,A,F9.4))') ' Hn(', nc4a,') = ', dn(nc4a), '  Hn(', nc4b,') = ', dn(nc4b)
  WRITE (UNIT=*, FMT='(A)') ' '

END SUBROUTINE Get_Classes

SUBROUTINE CnEc (nca, ncb)

  IMPLICIT NONE

  INTEGER (KIND=ki4), INTENT (IN) :: nca, ncb

  WRITE (0, '(A,2I5)') ' CnEc: nca, ncb = ', nca, ncb
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(SUM(ec(i,j,nca:ncb)),kr4)
    END DO
  END DO
WRITE (UNIT=50) buf

END SUBROUTINE CnEc

SUBROUTINE CnEp (nca, ncb)

  IMPLICIT NONE

  INTEGER (KIND=ki4), INTENT (IN) :: nca, ncb

  WRITE (0, '(A,2I5)') ' CnEp: nca, ncb = ', nca, ncb
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(SUM(ep(i,j,nca:ncb)),kr4)
    END DO
  END DO
WRITE (UNIT=50) buf

END SUBROUTINE CnEp

SUBROUTINE CnEt (nca, ncb)

  IMPLICIT NONE

  INTEGER (KIND=ki4), INTENT (IN) :: nca, ncb

  WRITE (0, '(A,2I5)') ' CnEt: nca, ncb = ', nca, ncb
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(SUM(et(i,j,nca:ncb)),kr4)
    END DO
  END DO
WRITE (UNIT=50) buf

END SUBROUTINE CnEt

SUBROUTINE CnPe (nca, ncb)

  IMPLICIT NONE

  INTEGER (KIND=ki4), INTENT (IN) :: nca, ncb

  WRITE (0, '(A,2I5)') ' CnPe: nca, ncb = ', nca, ncb
  DO j=1,Jmax
    DO i=1,Imax
      buf(i,j)=REAL(SUM(pe(i,j,nca:ncb)),kr4)
    END DO
  END DO
WRITE (UNIT=50) buf

END SUBROUTINE CnPe


SUBROUTINE Cm (nca, ncb)

  IMPLICIT NONE

  INTEGER (KIND=ki4), INTENT (IN) :: nca, ncb

  WRITE (0, '(A,2I5)') ' Cm: nca, ncb = ', nca, ncb
  !Class n
  DO k=Kmax,1,-1
    buf=0.0
    DO n=nca,ncb
      DO j=1,Jmax
        DO i=1,Imax
          etnn=0.5*cn2(n)*(un(i,j,n)*un(i,j,n)+  &
               vn(i,j,n)*vn(i,j,n)+hn(i,j,n)*hn(i,j,n))
          buf(i,j)=buf(i,j)+REAL(etnn*fn(k,n)*fn(k,n),kr4)
        END DO
      END DO
    END DO
    CALL Global_Mean (buf,bgm)
    WRITE (UNIT=70) buf
    WRITE (UNIT=75) bgm
  END DO

END SUBROUTINE Cm

SUBROUTINE CnCm (ncaa, ncab, ncba, ncbb)

  IMPLICIT NONE

  INTEGER (KIND=ki4), INTENT (IN) :: ncaa, ncab, ncba, ncbb

  WRITE (0, '(A,4I5)') ' CnCm: ncaa, ncab, ncba ncbb = ', ncaa, ncab, ncba, ncbb
  !Interaction Classes m with n
  DO k=Kmax,1,-1
    buf=0.0
    DO n=ncaa,ncab
      DO m=ncba,ncbb
        DO j=1,Jmax
          DO i=1,Imax
            etnm=gv*SQRT(dn(n)*dn(m))*(un(i,j,n)*un(i,j,m)+  &
                 vn(i,j,n)*vn(i,j,m)+hn(i,j,n)*hn(i,j,m))
            buf(i,j)=buf(i,j)+REAL(etnm*fn(k,n)*fn(k,m),kr4)
          END DO
        END DO
      END DO
    END DO
    CALL Global_Mean (buf,bgm)
    WRITE (UNIT=70) buf
    WRITE (UNIT=75) bgm
  END DO

END SUBROUTINE CnCm

END PROGRAM Vertical_Decomposition