PROGRAM cdfeke_KIEL
  !!======================================================================
  !!                 ***  PROGRAM cdfeke_KIEL   ***
  !!=====================================================================
  !!  ** Purpose : Compute Eddy Kinetic Energy 
  !!
  !!  ** Method  : Use 5d gridU gridU2, gridV gridV2 files
  !!               Velocities are interpolated both on T points
  !!               and the variance is computed. If -mke (-nie) option is used
  !!               the program also outputs MKE ("Near Inertial Energy") field
  !!
  !! History : pre  : 11/2004  : J.M. Molines : Original code
  !!           2.1  : 04/2005  : J.M. Molines : use modules
  !!           3.0  : 12/2010  : J.M. Molines : Doctor norm + Lic.
  !!           3.0.1: 11/2015  : J.K. Rieck   : GEOMAR-style EKE
  !!           3.0.2: 11/2015  : J.K. Rieck   : prescribe mean U/V
  !!----------------------------------------------------------------------
  USE cdfio 
  USE modcdfnames
  !!----------------------------------------------------------------------
  !! CDFTOOLS_3.0 , MEOM 2011
  !! $Id: cdfeke.f90 819 2015-05-23 07:57:41Z molines $
  !! Copyright (c) 2010, J.-M. Molines
  !! Software governed by the CeCILL licence (Licence/CDFTOOLSCeCILL.txt)
  !!----------------------------------------------------------------------
  IMPLICIT NONE

  INTEGER(KIND=4)                            :: ji, jj, jk, jt, jtm ! dummy loop index
  INTEGER(KIND=4)                            :: narg, iargc         ! command line browsing
  INTEGER(KIND=4)                            :: ijarg, ixtra, nfree ! command line browsing
  INTEGER(KIND=4)                            :: npiglo, npjglo      ! size of the domain (horiz)
  INTEGER(KIND=4)                            :: npk, npt, nt        ! size of the domain vert and time
  INTEGER(KIND=4)                            :: ncout, nptmean      ! ncid of output file
  INTEGER(KIND=4)                            :: ierr                ! Error status
  INTEGER(KIND=4)                            :: ivar                ! variable counter
  INTEGER(KIND=4)                            :: ip_eke, ip_mke      ! variable index
  INTEGER(KIND=4)                            :: ip_nie              ! variable index
  INTEGER(KIND=4), DIMENSION(3)              :: ipk, id_varout      ! 

  REAL(KIND=8)                               :: ua, va              ! working arrays
  REAL(KIND=4), DIMENSION(:),    ALLOCATABLE :: tim                 ! time variable
  REAL(KIND=8), DIMENSION (:,:), ALLOCATABLE :: uc, vc, u2, v2      ! velocities etc...
  REAL(KIND=8), DIMENSION (:,:), ALLOCATABLE :: eke, nie, um, vm    ! Eddy Kinetic Energy and Near Inertial Energy
  REAL(KIND=8), DIMENSION (:,:), ALLOCATABLE :: umean, vmean        !
  REAL(KIND=8), DIMENSION (:,:), ALLOCATABLE :: rmke                ! Mean Kinetic energy

  CHARACTER(LEN=256)                         :: cf_out='eke.nc'     ! file name
  CHARACTER(LEN=256)                         :: cf_ufil, cf_u2fil   ! file name
  CHARACTER(LEN=256)                         :: cf_vfil, cf_v2fil   !   "
  CHARACTER(LEN=256)                         :: cf_umfil, cf_vmfil  !   "
  CHARACTER(LEN=256)                         :: cf_tfil             !   "
  CHARACTER(LEN=256)                         :: cdum                ! dummy character variable

  TYPE(variable), DIMENSION(3)               :: stypvar             !

  LOGICAL                                    :: lchk                ! checking files existence
  LOGICAL                                    :: lperio=.FALSE.      ! checking E-W periodicity
  LOGICAL                                    :: leke=.TRUE.         ! compute EKE
  LOGICAL                                    :: lmke=.FALSE.        ! compute MKE
  LOGICAL                                    :: lnie=.FALSE.        ! commpute NIE
  LOGICAL                                    :: lnc4=.FALSE.        ! netcdf4 output (cunking and deflation)
  LOGICAL                                    :: omean=.FALSE.       ! calculate EKE with prescribed mean
  !!----------------------------------------------------------------------
  CALL ReadCdfNames()

  !!  Read command line
  narg= iargc()
  IF ( narg == 0 ) THEN
     PRINT *,' usage : cdfeke U-file [U2-file] [U-mean-file] V-file [V2-file] ...'
     PRINT *,'            ... [V-mean-file] T-file [-mke] [-nie] [-nc4] ...'
     PRINT *,'            ... [-omean] [-o output_file]'
     PRINT *,'      '
     PRINT *,'     PURPOSE :'
     PRINT *,'        Compute the Eddy Kinetic Energy from 5d mean values'
     PRINT *,'        and yearly mean of velocity components.'
     PRINT *,'        Preferably 5d data in files containing 1 year worth'
     PRINT *,'        of data should be used. EKE is calculated with'
     PRINT *,'        respect to U and V averaged over the length of'
     PRINT *,'        the time variable in the input files. All files'
     PRINT *,'        must share same time,depth,lon,lat dimensions!'
     PRINT *,'        (u-mean and v-mean can have a time dimension'
     PRINT *,'        of size 1, which is then used at every time step)' 
     PRINT *,'      '
     PRINT *,'     ARGUMENTS : both ''General Use'' or ''Reduced Use'' are acceptable'
     PRINT *,'       '
     PRINT *,'      * General Use: 3 files are given in argument, EKE is computed. '
     PRINT *,'       U-file  : gridU type file with 5d U component.' 
     PRINT *,'       V-file  : gridV type file with 5d V component.' 
     PRINT *,'       T-file  : any gridT or gridT2 (smaller) file, used for EKE header.'
     PRINT *,'       '
     PRINT *,'      * Extended Use: 5 files are given in argument, and EKE and is computed,'
     PRINT *,'       additionally "Near-Inertial Energy" NIE is computed.'
     PRINT *,'       U-file  : gridU type file with 5d U component.' 
     PRINT *,'       U2-file : gridU2 type file with 5d U2 component.' 
     PRINT *,'       V-file  : gridV type file with 5d V component.' 
     PRINT *,'       V2-file : gridV2 type file with 5d V2 component.' 
     PRINT *,'       T-file  : any gridT or gridT2 (smaller) file, used for EKE header.'
     PRINT *,'             '
     PRINT *,'      * Advanced Use: 5 or 7 files are given in argument, EKE is'
     PRINT *,'       calculated with respect to mean U and V given in U-mean and'
     PRINT *,'       V-mean files'
     PRINT *,'       !!!!!!!! Only works with -omean option !!!!!!!!!'
     PRINT *,'       U-file      : gridU type file with 5d U component'
     PRINT *,'       U-mean-file : gridU type file with mean U component'
     PRINT *,'       U2-file     : gridU2 type file with 5d U2 component. [optional]'
     PRINT *,'       V-file      : gridV type file with 5d V component'
     PRINT *,'       V-mean-file : gridV type file with mean V component'
     PRINT *,'       V2-file     : gridV2 type file with 5d V2 component. [optional]'
     PRINT *,'       T-file  : any gridT or gridT2 (smaller) file, used for EKE header.'
     PRINT *,'              '
     PRINT *,'     OPTION :'
     PRINT *,'       -mke  : output MKE field together with EKE. '
     PRINT *,'       -nie  : output NIE field together with EKE. '
     PRINT *,'       -omean  : compute EKE with respect to precribed mean. '
     PRINT *,'       -nc4  : allow netcdf4 output with compression and chunking.'
     PRINT *,'       -o output file : specify output file name instead of ', TRIM(cf_out)
     PRINT *,'      '
     PRINT *,'     REQUIRED FILES :'
     PRINT *,'        none'
     PRINT *,'      '
     PRINT *,'     OUTPUT : '
     PRINT *,'       netcdf file : ', TRIM(cf_out) , ' unless -o option in use.'
     PRINT *,'         variables : voeke (m2/s2)'
     PRINT *,'         variables : vomke (m2/s2) if required'
     PRINT *,'         variables : vonie (m2/s2) if required'
     STOP
  ENDIF
  !!
  !! Initialisation from 1st file (all file are assume to have the same geometry)
  cf_u2fil='none'
  cf_v2fil='none'
  cf_umfil='none'
  cf_vmfil='none'
  ip_eke=0
  ip_mke=0
  ip_nie=0
  ijarg = 1 ; ixtra = 0

  ! read arguments in 2 step (1) determine the number of 'free' arguments : 3 or 5 or 7
  !                          (2) according to the number of free arguments, set file names
  ! ( 1 ) 
  DO WHILE ( ijarg <= narg )
     CALL getarg(ijarg, cdum ) ; ijarg = ijarg + 1
     SELECT CASE ( cdum )
     CASE ( '-mke' )
         lmke = .true.
     CASE ( '-nie' )
         lnie = .true.
     CASE ( '-omean' )
         omean = .true.
     CASE ( '-nc4' )
         lnc4 = .true.
     CASE ( '-o'   )
         CALL getarg( ijarg, cf_out) ; ijarg = ijarg + 1
     CASE DEFAULT
       ! count xtra arguments ( ie those without key )
       ixtra = ixtra + 1
     END SELECT
  ENDDO
  nfree = ixtra
  ijarg = 1 ; ixtra = 0  ! reset ijarg for second step

  IF ( nfree /= 3 .AND. nfree /= 5 .AND. nfree /= 7) THEN
    PRINT *, ' +++ ERROR : not the correct number of free arguments'
    PRINT *, '            Likely a wrong option or missing file name'
    STOP
  ENDIF

  ! ( 2 )
  DO WHILE ( ijarg <= narg )
     CALL getarg(ijarg, cdum ) ; ijarg = ijarg + 1
     SELECT CASE ( cdum )
     CASE ( '-mke' )
     CASE ( '-nie' )
     CASE ( '-omean' )
     CASE ( '-nc4' )
     CASE ( '-o'   )
         ijarg = ijarg + 1
     CASE DEFAULT
       ixtra = ixtra + 1 
       SELECT CASE ( nfree )
       CASE ( 3 )
          IF ( lnie .OR. omean ) THEN
            PRINT *,'ERROR: -nie or -omean option with only 3 input files'
            PRINT *,'         is not possible!'
            STOP
          ENDIF
          IF ( ixtra == 1 ) cf_ufil = cdum
          IF ( ixtra == 2 ) cf_vfil = cdum
          IF ( ixtra == 3 ) cf_tfil = cdum
          !lmke = .true.
          leke = .true.
          lnie = .false.
          omean = .false.
       CASE ( 5 )
          IF ( omean ) THEN
            IF ( lnie ) THEN
              PRINT *,'ERROR: -nie AND -omean option with only 5 input files'
              PRINT *,'         is not possible!'
              STOP
            ENDIF
            IF ( ixtra == 1 ) cf_ufil  = cdum
            IF ( ixtra == 2 ) cf_umfil = cdum
            IF ( ixtra == 3 ) cf_vfil  = cdum
            IF ( ixtra == 4 ) cf_vmfil = cdum
            IF ( ixtra == 5 ) cf_tfil  = cdum
            lnie = .false.
            !lmke = .true.
            leke = .true.
            !omean = .true.
          ELSE
            IF ( ixtra == 1 ) cf_ufil  = cdum
            IF ( ixtra == 2 ) cf_u2fil = cdum
            IF ( ixtra == 3 ) cf_vfil  = cdum
            IF ( ixtra == 4 ) cf_v2fil = cdum
            IF ( ixtra == 5 ) cf_tfil  = cdum
            !lnie = .true.
            !lmke = .true.
            !leke = .true.
            omean = .false.
          ENDIF
        CASE ( 7 )
            IF ( ixtra == 1 ) cf_ufil  = cdum
            IF ( ixtra == 2 ) cf_u2fil = cdum
            IF ( ixtra == 3 ) cf_umfil = cdum
            IF ( ixtra == 4 ) cf_vfil  = cdum
            IF ( ixtra == 5 ) cf_v2fil = cdum
            IF ( ixtra == 6 ) cf_vmfil = cdum
            IF ( ixtra == 7 ) cf_tfil  = cdum
            lnie = .true.
            lmke = .true.
            leke = .true.
            omean = .true.
       END SELECT
     END SELECT
  ENDDO

  lchk =           chkfile (cf_ufil )
  lchk = lchk .OR. chkfile (cf_vfil )
  lchk = lchk .OR. chkfile (cf_tfil )
  IF (  lnie ) THEN
    lchk = lchk .OR. chkfile (cf_u2fil)
    lchk = lchk .OR. chkfile (cf_v2fil)
  ENDIF
  IF ( omean ) THEN
    lchk = lchk .OR. chkfile (cf_umfil)
    lchk = lchk .OR. chkfile (cf_vmfil)
  ENDIF

  IF ( lchk ) STOP ! missing files

  npiglo = getdim (cf_ufil,cn_x)
  npjglo = getdim (cf_ufil,cn_y)
  npk    = getdim (cf_ufil,cn_z)
  npt    = getdim (cf_ufil,cn_t)
  nptmean = getdim (cf_umfil, cn_t)

  IF ( npk == 0 ) npk=1 ! assume 1 level at least

  stypvar(1)%ichunk = (/ npiglo, npjglo, 1, 1 /) 
  stypvar(2)%ichunk = (/ npiglo, npjglo, 1, 1 /) 
  stypvar(3)%ichunk = (/ npiglo, npjglo, 1, 1 /) 
  
  ivar = 1
  IF ( leke ) THEN 
    ip_eke=ivar
    ipk(ip_eke)                       = npk
    stypvar(ip_eke)%cname             = 'voeke'
    stypvar(ip_eke)%cunits            = 'm2/s2'
    stypvar(ip_eke)%rmissing_value    = 0.
    stypvar(ip_eke)%valid_min         = 0.
    stypvar(ip_eke)%valid_max         = 10000.
    stypvar(ip_eke)%clong_name        = 'Eddy_Kinetic_Energy'
    stypvar(ip_eke)%cshort_name       = 'voeke'
    stypvar(ip_eke)%conline_operation = 'N/A'
    stypvar(ip_eke)%caxis             = 'TZYX'
    ivar = ivar + 1
  ENDIF

  IF ( lmke ) THEN
    ip_mke=ivar
    ipk(ip_mke)                       = npk
    stypvar(ip_mke)%cname             = 'vomke'
    stypvar(ip_mke)%cunits            = 'm2/s2'
    stypvar(ip_mke)%rmissing_value    = 0.
    stypvar(ip_mke)%valid_min         = 0.
    stypvar(ip_mke)%valid_max         = 10000.
    stypvar(ip_mke)%clong_name        = 'Mean_Kinetic_Energy'
    stypvar(ip_mke)%cshort_name       = 'vomke'
    stypvar(ip_mke)%conline_operation = 'N/A'
    stypvar(ip_mke)%caxis             = 'TZYX'
    ivar = ivar + 1
  ENDIF
  
  IF ( lnie ) THEN
    ip_nie=ivar
    ipk(ip_nie)                       = npk
    stypvar(ip_nie)%cname             = 'vonie'
    stypvar(ip_nie)%cunits            = 'm2/s2'
    stypvar(ip_nie)%rmissing_value    = 0.
    stypvar(ip_nie)%valid_min         = 0.
    stypvar(ip_nie)%valid_max         = 10000.
    stypvar(ip_nie)%clong_name        = 'Near_Inertial_Energy'
    stypvar(ip_nie)%cshort_name       = 'vonie'
    stypvar(ip_nie)%conline_operation = 'N/A'
    stypvar(ip_nie)%caxis             = 'TZYX'
  ENDIF

  ivar = MAX( ip_mke, ip_eke, ip_nie) ! set ivar to the effective number of variables to be output

  PRINT *, 'npiglo = ', npiglo
  PRINT *, 'npjglo = ', npjglo
  PRINT *, 'npk    = ', npk
  PRINT *, 'npt    = ', npt

  ALLOCATE( uc(npiglo,npjglo), vc(npiglo,npjglo) )
  ALLOCATE( tim(npt) )

  IF ( leke ) THEN
    ALLOCATE( eke(npiglo,npjglo)  )
    ALLOCATE( um(npiglo,npjglo), vm(npiglo,npjglo) )
    ALLOCATE( umean(npiglo,npjglo), vmean(npiglo,npjglo) )
    eke(:,:) = 0.e0
    umean(:,:) = 0.e0
    vmean(:,:) = 0.e0
  ENDIF
  IF ( lnie ) THEN
    ALLOCATE( nie(npiglo,npjglo)  )
    ALLOCATE( u2(npiglo,npjglo), v2(npiglo,npjglo) )
    nie(:,:) = 0.e0
  ENDIF
  IF (lmke ) THEN
    ALLOCATE( rmke(npiglo,npjglo)  )
    rmke(:,:) = 0.e0
  ENDIF

  ncout = create      (cf_out, cf_tfil, npiglo, npjglo, npk       , ld_nc4=lnc4 )
  ierr  = createvar   (ncout,  stypvar, ivar,   ipk,    id_varout , ld_nc4=lnc4 )
  ierr  = putheadervar(ncout,  cf_tfil, npiglo, npjglo, npk                     )

  ! check for E_W periodicity
  uc(:,:) = getvar(cf_tfil, cn_vlon2d, 1, npiglo, npjglo, ktime=1 )
  IF ( uc(1,1) ==  uc(npiglo-1,1) ) THEN 
     lperio = .TRUE. 
     PRINT *,' E-W periodicity detected '
  ENDIF

  DO jk = 1, npk
    IF ( omean ) THEN
    ELSE
    um(:,:) = 0.e0
    vm(:,:) = 0.e0
    nt = 0
    DO jtm = 1, npt
      nt = nt + 1
      uc(:,:) = getvar(cf_ufil,  cn_vozocrtx,               jk, npiglo, npjglo, ktime=jtm )
      vc(:,:) = getvar(cf_vfil,  cn_vomecrty,               jk, npiglo, npjglo, ktime=jtm )
      um(:,:) = um(:,:) + uc(:,:)*1.d0
      vm(:,:) = vm(:,:) + vc(:,:)*1.d0
    END DO
    umean(:,:) = um(:,:)/nt
    vmean(:,:) = vm(:,:)/nt
    ENDIF
    DO jt = 1, npt
      IF ( omean ) THEN
        IF ( nptmean == 1 ) THEN
           umean(:,:) = getvar(cf_umfil,cn_vozocrtx,        jk, npiglo, npjglo, ktime=1)
           vmean(:,:) = getvar(cf_vmfil,cn_vomecrty,        jk, npiglo, npjglo, ktime=1)
        ELSE
           umean(:,:) = getvar(cf_umfil,cn_vozocrtx,           jk, npiglo, npjglo, ktime=jt)
           vmean(:,:) = getvar(cf_vmfil,cn_vomecrty,           jk, npiglo, npjglo, ktime=jt)
        ENDIF
      ENDIF
      uc(:,:) = getvar(cf_ufil,  cn_vozocrtx,               jk, npiglo, npjglo, ktime=jt )
      vc(:,:) = getvar(cf_vfil,  cn_vomecrty,               jk, npiglo, npjglo, ktime=jt )
      IF ( lnie ) THEN
        u2(:,:) = getvar(cf_u2fil, TRIM(cn_vozocrtx)//'_sqd', jk ,npiglo, npjglo, ktime=jt )
        v2(:,:) = getvar(cf_v2fil, TRIM(cn_vomecrty)//'_sqd', jk ,npiglo, npjglo, ktime=jt )
      ENDIF

      ua = 0.e0 ; va = 0.e0 ; eke(:,:) = 0.e0
      IF ( leke ) THEN
        DO ji=2, npiglo
          DO jj=2,npjglo
              ua = 0.5* (((uc(ji,jj)-umean(ji,jj))*(uc(ji,jj)-umean(ji,jj))) + &
              &          ((uc(ji-1,jj)-umean(ji-1,jj))*(uc(ji-1,jj)-umean(ji-1,jj))))
              va = 0.5* (((vc(ji,jj)-vmean(ji,jj))*(vc(ji,jj)-vmean(ji,jj))) + &
              &          ((vc(ji,jj-1)-vmean(ji,jj-1))*(vc(ji,jj-1)-vmean(ji,jj-1))))
              eke(ji,jj) = 0.5 * ( ua + va )
          END DO
        END DO
      ENDIF
      
      ua = 0.e0 ; va = 0.e0 ; nie(:,:) = 0.e0
      IF ( lnie ) THEN
        DO ji=2, npiglo
          DO jj=2,npjglo
              ua = 0.5* ((u2(ji,jj)-uc(ji,jj)*uc(ji,jj))+ (u2(ji-1,jj)-uc(ji-1,jj)*uc(ji-1,jj)))
              va = 0.5* ((v2(ji,jj)-vc(ji,jj)*vc(ji,jj))+ (v2(ji,jj-1)-vc(ji,jj-1)*vc(ji,jj-1)))
              nie(ji,jj) = 0.5 * ( ua + va )
          END DO
        END DO
      ENDIF
      
      rmke(:,:) = 0.e0
      IF ( lmke ) THEN
        DO ji=2, npiglo
          DO jj=2,npjglo
              rmke(ji,jj)=  0.5* (0.5*( umean(ji,jj)*umean(ji,jj) + umean(ji-1,jj)*umean(ji-1,jj)) + &
              &                   0.5*( vmean(ji,jj)*vmean(ji,jj) + vmean(ji,jj-1)*vmean(ji,jj-1)) )
          END DO
        END DO
      ENDIF

      IF ( lperio ) eke(1,:) = eke(npiglo-1,:)
      IF ( leke ) THEN 
         IF ( lperio ) eke(1,:) = eke(npiglo-1,:)
         ierr=putvar(ncout,id_varout(ip_eke), eke,  jk ,npiglo, npjglo, ktime=jt )
      ENDIF
      IF ( lmke ) THEN 
         IF ( lperio ) rmke(1,:) = rmke(npiglo-1,:)
         ierr=putvar(ncout,id_varout(ip_mke), rmke, jk ,npiglo, npjglo, ktime=jt )
      ENDIF
      IF ( lnie ) THEN 
         IF ( lperio ) nie(1,:) = nie(npiglo-1,:)
         ierr=putvar(ncout,id_varout(ip_nie), nie, jk ,npiglo, npjglo, ktime=jt )
      ENDIF
    END DO
  END DO ! time loop

  tim  = getvar1d(cf_ufil, cn_vtimec, npt     )
  ierr = putvar1d(ncout,   tim,       npt, 'T')

  ierr = closeout(ncout)

END PROGRAM cdfeke_KIEL