Subversion Repositories lagranto.icon

! $RCSfile: utilities.f90,v $

! $Revision: 4.11 $ $Date: 2009/11/30 14:29:09 $

!+ Source module for utility routines

!==============================================================================

MODULE  utilities

!==============================================================================

!

! Description:

!   This module provides service utilities for the model. All routines are 

!   written in a manner that also other models can use it. That means:

!     - no routine uses other modules, except the declarations for the 

!       KIND-type parameter; the data access is by parameter list only

!     - no routine allocates dynamic memory; work space needed is

!       provided via the parameter list

!     - no derived data types are used

!

!   Routines (module procedures) currently contained:

!

!     - convert_month:

!       Converts a 3-character string abbreviation of a month into the number

!       of the month or vice versa.

!

!     - dfilt4:

!       Digital filter of length 4

!

!     - dfilt8:

!       Digital filter of length 8

!

!     - dolph:

!       Calculates the Dolph-Chebyshev window for the initialization

!

!     - elapsed_time:

!       Returns the elapsed wall-clock time in seconds since the last call.

!       On the first call the variables are only initialized. If no system

!       clock is present, an error-value will be returned

!

!     - get_utc_date:

!       Calculates the actual date using the date of the forecast-start and 

!       the number of timesteps performed.

!

!     - horizontal_filtering

!       horizontal filtering (at the moment especially for the pressure deviation)

!

!     - phirot2phi:

!       Converts phi from the rotated system to phi in the real

!       geographical system.

!

!     - phi2phirot:

!       Converts phi from the real geographical system to phi

!       in the rotated system.

!

!     - rlarot2rla:

!       Converts lambda from the rotated system to lambda in the real

!       geographical system.

!

!     - rla2rlarot:

!       Converts lambda from the real geographical system to lambda 

!       in the rotated system.

!

!     - sleve_split_oro

!       Decomposes a given topography field in a large-scale and a small-scale

!       part according to the definition of the SLEVE coordinate

!

!     - smoother:

!       Smoothes a 2-D field by applying digital filters

!

!     - tautsp:

!       Computes tension splines

!

!     - tautsp2D:

!       Computes tension splines for several columns

!

!     - to_upper:

!       Converts alphabetic characters from lower to upper case.

!

!     - uvrot2uv:

!       Converts the wind components u and v from the rotated system

!       to the real geographical system.

!

!     - uvrot2uv_vec:

!       the same as above, but for a whole 2D field (in vectorized form).

!

!     - uv2uvrot:

!       Converts the wind components u and v from the real geographical

!       system to the rotated system.

!

!     - uv2uvrot_vec:

!       the same as above, but for a whole 2D field (in vectorized form).

!

!     - uv2df:

!       Converts the wind components u and v to wind direction and speed.

!

!     - uv2df_vec:

!       the same as above, but for a whole 2D field (in vectorized form).

!

!

! Current Code Owner: DWD, Ulrich Schaettler

!  phone:  +49  69  8062 2739

!  fax:    +49  69  8062 3721

!  email:  ulrich.schaettler@dwd.de

!

! History:

! Version    Date       Name

! ---------- ---------- ----

! 1.1        1998/03/11 Ulrich Schaettler

!  Initial release

! 1.2        1998/03/30 Ulrich Schaettler

!  Introduction of subroutine dolph used during the initialization

! 1.9        1998/09/16 Guenther Doms

!  Introduction of a smoothing routine 'smoother' which uses digital

!  filters 'dfilt4' and 'dfilt8'.

! 1.10       1998/09/29 Ulrich Schaettler

!  Routine remark eliminated and put to parallel_utilities.

!  Routines uv2uvrot and uv2df introduced

! 1.16       1998/11/02 Guenther Doms

!  Correction of filter processing in routine 'smoother'.

! 1.29       1999/05/11 Ulrich Schaettler

!  Adaptations to use this module also in GME2LM

! 1.32       1999/08/24 Guenther Doms

!  some _ireals declarations added.

! 2.8        2001/07/06 Ulrich Schaettler

!  Added new subroutines tautsp2D, uv2uvrot_vec and uvrot2uv_vec for 

!  vectorization

! 2.14       2002/02/15 Ulrich Schaettler

!  Correction and adaptations in tautsp2D (analogous to GME2LM)

!  Added new subroutine dc_topo for the SLEVE coordinate

! 2.17       2002/05/08 Ulrich Schaettler

!  Modifications for performing the filtering in irealgrib-format

! 2.18       2002/07/16 Guenther Doms

!  Corrections for the rotation of the wind components from or to the

!  geographical coordinate system.

! 3.3        2003/04/22 Christoph Schraff

!  Introduction of subroutines 'convert_month' and 'to_upper' (for GPS data).

! 3.6        2003/12/11 Ulrich Schaettler

!  Eliminated Subroutine istringlen (use F90 intrinsic LEN_TRIM instead)

! 3.13       2004/12/03 Ulrich Schaettler

!  Eliminated dependency on data_io (put irealgrib to data_parameters)

!  New SR horizontal_filtering (from INT2LM);

!  Renamed SR dc_topo to sleve_split_oro

! 3.14       2005/01/25 Ulrich Schaettler

!  New filter routine smooth9 for new type of Rayleigh damping (Lucio Torrisi)

!  Changes in horizontal_filtering (Jochen Foerstner)

! 3.15       2005/03/03 Ulrich Schaettler

!  Replaced FLOAT by REAL

! 3.16       2005/07/22 Ulrich Schaettler

!  Bug correction in the call to intrinsic function REAL

! 3.18       2006/03/03 Ulrich Schaettler

!  Introduced idouble/isingle as KIND parameters instead of ireals/irealgrib

!  in the generic formulation of some routines (dfilt4, dfilt8, smoother)

!  Changed get_utc_date to include also a climatological year with 360 days

! 3.21       2006/12/04 Burkhardt Rockel, Lucio Torrisi, Jochen Foerstner

!  Added polgam in transformation function rla2rlarot

!  polgam is not used as optional parameter any more

!  Some adaptations in smooth9 for itype_spubc=2

!  Some modifications in horizontal_filtering

!  Function uv2df_vec introduced. (C. Schraff)

! V3_23        2007/03/30 Ulrich Schaettler

!  Declared some constant variables as parameters to allow inlining on

!  some platforms

!  Changed computation of acthour in get_utc_date

! V3_24        2007/04/26 Ulrich Schaettler

!  Bug correction in computation of acthour in get_utc_date

! V4_1         2007/12/04 Ulrich Schaettler

!  Introduced parameter myid to sleve_split_oro (is called from all PEs in INT2LM)

! V4_4         2008/07/16 Ulrich Schaettler

!  Adapted a debug printout in SR tautsp2D

!  Changed NL parameter lyear_360 to itype_calendar, to have several options

!  Vectorized SR horizontal_filtering by changing some loops

!  Treatment of very small values for spline interpolation in tautsp2D

! V4_8         2009/02/16 Ulrich Schaettler (Andy Dobler)

!  Corrected leap year calculation for centuries in Gregorian calendar

! @VERSION@    @DATE@     <Your name>

!  <Modification comments>         

!

! Code Description:

! Language: Fortran 90.

! Software Standards: "European Standards for Writing and

! Documenting Exchangeable Fortran 90 Code".

!==============================================================================

!

! Declarations:

!

! Modules used:

USE data_parameters , ONLY :   &

    ireals,    & ! KIND-type parameter for real variables

    iintegers, & ! KIND-type parameter for standard integer variables

    irealgrib, & ! KIND-type parameter for real variables in the grib library

    idouble,   & ! KIND-type parameter for double precision real variables

    isingle      ! KIND-type parameter for single precision real variables

!==============================================================================

IMPLICIT NONE

!==============================================================================

! Interface Blocks

INTERFACE smoother

  MODULE PROCEDURE                        &

    smoother_double,                      &

    smoother_single

END INTERFACE

INTERFACE dfilt4

  MODULE PROCEDURE                        &

    dfilt4_double,                        &

    dfilt4_single

END INTERFACE

INTERFACE dfilt8

  MODULE PROCEDURE                        &

    dfilt8_double,                        &

    dfilt8_single

END INTERFACE

!==============================================================================

CONTAINS

!==============================================================================

SUBROUTINE convert_month ( MonthString, MonthNumber, ind )

!-------------------------------------------------------------------------------

!

! Description:

!   Convert 3-chr Month string to number (ind > 0) or vice-versa (ind <= 0).

!     If ind > 0 and input string not valid, MonthNumber will be 0.

!     If ind <=0 and input month number invalid, MonthString will be 'XXX'.

! Method:

!     Uses subroutine 'to_upper'.

!-------------------------------------------------------------------------------

  IMPLICIT NONE

! Subroutine arguments:

! --------------------

  CHARACTER (LEN=3)       , INTENT(INOUT) :: MonthString  ! 3-chr Month name

  INTEGER (KIND=iintegers), INTENT(INOUT) :: MonthNumber  ! Month number (1-12)

  INTEGER (KIND=iintegers), INTENT(IN)    :: ind          ! > 0 : chr --> no.

                                                          ! <=0 : no  --> chr

! Local parameters:

! ----------------

  CHARACTER (LEN=36), PARAMETER ::  &

    MonthNames = "JANFEBMARAPRMAYJUNJULAUGSEPOCTNOVDEC"

! Local variables

! ---------------

  CHARACTER (LEN=3)             :: Month

  INTEGER (KIND=iintegers)      :: idx

!

!------------ End of header ----------------------------------------------------

  IF ( ind > 0 ) THEN

! ----- String to number -----

    Month = MonthString

    CALL to_upper ( Month )

    idx = INDEX ( MonthNames, Month )

    IF ( MOD ( idx-1, 3 ) == 0 ) THEN

      MonthNumber = ( idx + 2 ) / 3

    ELSE

      MonthNumber = 0

    END IF

  ELSE

! ----- Number to string -----

    IF ( MonthNumber >= 1 .AND. &

         MonthNumber <= 12 ) THEN

      idx = MonthNumber * 3 - 2

      MonthString = MonthNames(idx:idx+2)

    ELSE

      MonthString = "XXX"

    END IF

  END IF

END SUBROUTINE convert_month

!------------------------------------------------------------------------------

!==============================================================================

!==============================================================================

!+ Defines all subroutines for the generic routine dfilt4

!------------------------------------------------------------------------------

!

!  SUBROUTINE dfilt4 (fin, idim, fhelp, fout, nfilt)

!

!------------------------------------------------------------------------------

!

! Description:

!   This routine smoothes an arbitrary field (fin) of length idim by applying

!   a digital filters of length nlength 4 nfilt times. The filterd field

!   is written on fout.

!

! Method:

!   Digital filter according to Shapiro

!

!------------------------------------------------------------------------------

!+ Implementation for double precision

!------------------------------------------------------------------------------

SUBROUTINE dfilt4_double (fin, idim, fhelp, fout, nfilt)

!------------------------------------------------------------------------------

! Parameter list:

INTEGER (KIND=iintegers), INTENT (IN)          ::    &

  idim,           & ! Dimension of the field

  nfilt             ! Number of iterative filerings

REAL (KIND=idouble),   INTENT (IN)          ::    &

  fin (idim)        ! input field (unfilterd)

REAL (KIND=idouble),   INTENT (OUT)         ::    &

  fout (idim)       ! smoothed output field (filtered)

REAL (KIND=idouble),   INTENT (INOUT)       ::    &

  fhelp(idim)       ! additional storage supplied by the calling routine

! Local variables

INTEGER (KIND=iintegers) ::    &

  i,m,            & ! loop indicees

  nf_o2             ! nfilt/2

REAL (KIND=idouble)      ::    & 

  fw(5)             ! filter weights

!------------------------------------------------------------------------------

  DATA fw / -0.00390625_idouble, 0.03125_idouble, -0.109375_idouble,     &

             0.21875_idouble,    0.7265625_idouble /

! begin subroutine dfilt4_double

  nf_o2 = (nfilt+1)/2

  fout (:) = fin(:)

  fhelp(:) = fin(:)

  DO i = 2, idim-1

    fhelp(i) =   0.15_idouble*fout (i-1) + 0.7_idouble*fout (i)    &

               + 0.15_idouble*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.15_idouble*fhelp(i-1) + 0.7_idouble*fhelp(i)    &

               + 0.15_idouble*fhelp(i+1)

  ENDDO

  DO m = 1, nf_o2

    DO i = 5, idim-4

      fhelp(i) =  fw(5)*fout(i) &

                + fw(4)*(fout(i-1)+fout(i+1)) + fw(3)*(fout(i-2)+fout(i+2)) &

                + fw(2)*(fout(i-3)+fout(i+3)) + fw(1)*(fout(i-4)+fout(i+4))  

    ENDDO

    DO i = 5, idim-4

      fout(i) = fw(5)*fhelp(i) &

              + fw(4)*(fhelp(i-1)+fhelp(i+1)) + fw(3)*(fhelp(i-2)+fhelp(i+2)) &

              + fw(2)*(fhelp(i-3)+fhelp(i+3)) + fw(1)*(fhelp(i-4)+fhelp(i+4))  

    ENDDO

  ENDDO

  DO i = 2, idim-1

    fhelp(i) =   0.15_idouble*fout (i-1) + 0.7_idouble*fout (i)   &

               + 0.15_idouble*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.15_idouble*fhelp(i-1) + 0.7_idouble*fhelp(i)   &

               + 0.15_idouble*fhelp(i+1)

  ENDDO

END SUBROUTINE dfilt4_double

!------------------------------------------------------------------------------

!+ Implementation for single precision

!------------------------------------------------------------------------------

SUBROUTINE dfilt4_single (fin, idim, fhelp, fout, nfilt)

!------------------------------------------------------------------------------

! Parameter list:

INTEGER (KIND=iintegers), INTENT (IN)          ::    &

  idim,           & ! Dimension of the field

  nfilt             ! Number of iterative filerings

REAL (KIND=isingle),   INTENT (IN)          ::    &

  fin (idim)        ! input field (unfilterd)

REAL (KIND=isingle),   INTENT (OUT)         ::    &

  fout (idim)       ! smoothed output field (filtered)

REAL (KIND=isingle),   INTENT (INOUT)       ::    &

  fhelp(idim)       ! additional storage supplied by the calling routine

! Local variables

INTEGER (KIND=iintegers) ::    &

  i,m,            & ! loop indicees

  nf_o2             ! nfilt/2

REAL (KIND=isingle)      ::    & 

  fw(5)             ! filter weights

!------------------------------------------------------------------------------

  DATA fw / -0.00390625_isingle, 0.03125_isingle, -0.109375_isingle,     &

             0.21875_isingle,    0.7265625_isingle /

! begin subroutine dfilt4_single

  nf_o2 = (nfilt+1)/2

  fout (:) = fin(:)

  fhelp(:) = fin(:)

  DO i = 2, idim-1

    fhelp(i) =   0.15_isingle*fout (i-1) + 0.7_isingle*fout (i)    &

               + 0.15_isingle*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.15_isingle*fhelp(i-1) + 0.7_isingle*fhelp(i)    &

               + 0.15_isingle*fhelp(i+1)

  ENDDO

  DO m = 1, nf_o2

    DO i = 5, idim-4

      fhelp(i) =  fw(5)*fout(i) &

                + fw(4)*(fout(i-1)+fout(i+1)) + fw(3)*(fout(i-2)+fout(i+2)) &

                + fw(2)*(fout(i-3)+fout(i+3)) + fw(1)*(fout(i-4)+fout(i+4))  

    ENDDO

    DO i = 5, idim-4

      fout(i) = fw(5)*fhelp(i) &

              + fw(4)*(fhelp(i-1)+fhelp(i+1)) + fw(3)*(fhelp(i-2)+fhelp(i+2)) &

              + fw(2)*(fhelp(i-3)+fhelp(i+3)) + fw(1)*(fhelp(i-4)+fhelp(i+4))  

    ENDDO

  ENDDO

  DO i = 2, idim-1

    fhelp(i) =   0.15_isingle*fout (i-1) + 0.7_isingle*fout (i)   &

               + 0.15_isingle*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.15_isingle*fhelp(i-1) + 0.7_isingle*fhelp(i)   &

               + 0.15_isingle*fhelp(i+1)

  ENDDO

END SUBROUTINE dfilt4_single

!------------------------------------------------------------------------------

!==============================================================================

!==============================================================================

!+ Defines all subroutines for the generic routine dfilt8

!------------------------------------------------------------------------------

!

! SUBROUTINE dfilt8 (fin, idim, fhelp, fout, nfilt)

!

!------------------------------------------------------------------------------

!

! Description:

!   This routine smoothes an arbitrary field (fin) of length idim by applying

!   a digital filters of length nlength 8 nfilt times. The filterd field

!   is written on fout.

!

! Method:

!   Digital filter according to Shapiro

!

!------------------------------------------------------------------------------

!+ Implementation for double precision

!------------------------------------------------------------------------------

SUBROUTINE dfilt8_double (fin, idim, fhelp, fout, nfilt)

!------------------------------------------------------------------------------

! Parameter list:

INTEGER (KIND=iintegers), INTENT (IN)          ::    &

  idim,           & ! Dimension of the field

  nfilt             ! Number of iterative filerings

REAL (KIND=idouble),   INTENT (IN)          ::    &

  fin (idim)        ! input field (unfilterd)

REAL (KIND=idouble),   INTENT (OUT)         ::    &

  fout (idim)       ! smoothed output field (filtered)

REAL (KIND=idouble),   INTENT (INOUT)       ::    &

  fhelp(idim)       ! additional storage supplied by the calling routine

! Local variables

INTEGER (KIND=iintegers) ::    &

  i,m,            & ! loop indicees

  nf_o2             ! nfilt/2

REAL (KIND=idouble)         ::  & 

  fw(9)  ! filter weights

!------------------------------------------------------------------------------

DATA fw /-0.0000152590_idouble,  0.0002441406_idouble, -0.0018310546_idouble, &

          0.0085449218_idouble, -0.0277709960_idouble,  0.0666503906_idouble, &

         -0.1221923828_idouble,  0.1745605469_idouble,  0.8036193848_idouble /

! begin subroutine dfilt8_double

  nf_o2 = (nfilt+1)/2

  fout (:) = fin(:)

  fhelp(:) = fin(:)

  DO i = 2, idim-1

    fhelp(i) =   0.25_idouble*fout (i-1) + 0.5_idouble*fout (i)     &

               + 0.25_idouble*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.25_idouble*fhelp(i-1) + 0.5_idouble*fhelp(i)     &

               + 0.25_idouble*fhelp(i+1)

  ENDDO

  DO m = 1, nf_o2

    DO i = 9, idim-8

      fhelp(i) = fw(9)*fout(i) &

               + fw(8)*(fout(i-1)+fout(i+1)) + fw(7)*(fout(i-2)+fout(i+2)) &

               + fw(6)*(fout(i-3)+fout(i+3)) + fw(5)*(fout(i-4)+fout(i+4)) &

               + fw(4)*(fout(i-5)+fout(i+5)) + fw(3)*(fout(i-6)+fout(i+6)) &

               + fw(2)*(fout(i-7)+fout(i+7)) + fw(1)*(fout(i-8)+fout(i+8))  

    ENDDO

    DO i = 9, idim-8

      fout(i) = fw(9)*fhelp(i) &

              + fw(8)*(fhelp(i-1)+fhelp(i+1)) + fw(7)*(fhelp(i-2)+fhelp(i+2)) &

              + fw(6)*(fhelp(i-3)+fhelp(i+3)) + fw(5)*(fhelp(i-4)+fhelp(i+4)) &

              + fw(4)*(fhelp(i-5)+fhelp(i+5)) + fw(3)*(fhelp(i-6)+fhelp(i+6)) &

              + fw(2)*(fhelp(i-7)+fhelp(i+7)) + fw(1)*(fhelp(i-8)+fhelp(i+8))

    ENDDO

  ENDDO

  DO i = 2, idim-1

    fhelp(i) =   0.25_idouble*fout (i-1) + 0.5_idouble*fout (i)    &

               + 0.25_idouble*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.25_idouble*fhelp(i-1) + 0.5_idouble*fhelp(i)    &

               + 0.25_idouble*fhelp(i+1)

  ENDDO

END SUBROUTINE dfilt8_double

!------------------------------------------------------------------------------

!+ Implementation for single precision

!------------------------------------------------------------------------------

SUBROUTINE dfilt8_single (fin, idim, fhelp, fout, nfilt)

!------------------------------------------------------------------------------

! Parameter list:

INTEGER (KIND=iintegers), INTENT (IN)          ::    &

  idim,           & ! Dimension of the field

  nfilt             ! Number of iterative filerings

REAL (KIND=isingle),   INTENT (IN)          ::    &

  fin (idim)        ! input field (unfilterd)

REAL (KIND=isingle),   INTENT (OUT)         ::    &

  fout (idim)       ! smoothed output field (filtered)

REAL (KIND=isingle),   INTENT (INOUT)       ::    &

  fhelp(idim)       ! additional storage supplied by the calling routine

! Local variables

INTEGER (KIND=iintegers) ::    &

  i,m,            & ! loop indicees

  nf_o2             ! nfilt/2

REAL (KIND=isingle)         ::  & 

  fw(9)  ! filter weights

!------------------------------------------------------------------------------

DATA fw /-0.0000152590_isingle,  0.0002441406_isingle, -0.0018310546_isingle, &

          0.0085449218_isingle, -0.0277709960_isingle,  0.0666503906_isingle, &

         -0.1221923828_isingle,  0.1745605469_isingle,  0.8036193848_isingle /

! begin subroutine dfilt8_single

  nf_o2 = (nfilt+1)/2

  fout (:) = fin(:)

  fhelp(:) = fin(:)

  DO i = 2, idim-1

    fhelp(i) =   0.25_isingle*fout (i-1) + 0.5_isingle*fout (i)     &

               + 0.25_isingle*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.25_isingle*fhelp(i-1) + 0.5_isingle*fhelp(i)     &

               + 0.25_isingle*fhelp(i+1)

  ENDDO

  DO m = 1, nf_o2

    DO i = 9, idim-8

      fhelp(i) = fw(9)*fout(i) &

               + fw(8)*(fout(i-1)+fout(i+1)) + fw(7)*(fout(i-2)+fout(i+2)) &

               + fw(6)*(fout(i-3)+fout(i+3)) + fw(5)*(fout(i-4)+fout(i+4)) &

               + fw(4)*(fout(i-5)+fout(i+5)) + fw(3)*(fout(i-6)+fout(i+6)) &

               + fw(2)*(fout(i-7)+fout(i+7)) + fw(1)*(fout(i-8)+fout(i+8))  

    ENDDO

    DO i = 9, idim-8

      fout(i) = fw(9)*fhelp(i) &

              + fw(8)*(fhelp(i-1)+fhelp(i+1)) + fw(7)*(fhelp(i-2)+fhelp(i+2)) &

              + fw(6)*(fhelp(i-3)+fhelp(i+3)) + fw(5)*(fhelp(i-4)+fhelp(i+4)) &

              + fw(4)*(fhelp(i-5)+fhelp(i+5)) + fw(3)*(fhelp(i-6)+fhelp(i+6)) &

              + fw(2)*(fhelp(i-7)+fhelp(i+7)) + fw(1)*(fhelp(i-8)+fhelp(i+8))

    ENDDO

  ENDDO

  DO i = 2, idim-1

    fhelp(i) =   0.25_isingle*fout (i-1) + 0.5_isingle*fout (i)    &

               + 0.25_isingle*fout (i+1)

  ENDDO

  DO i = 2, idim-1

    fout (i) =   0.25_isingle*fhelp(i-1) + 0.5_isingle*fhelp(i)    &

               + 0.25_isingle*fhelp(i+1)

  ENDDO

END SUBROUTINE dfilt8_single

!==============================================================================

!==============================================================================

!------------------------------------------------------------------------------

SUBROUTINE dolph (deltat, taus, m, window, t, time, time2, w, w2)

!------------------------------------------------------------------------------

!

! Description:

!  Calculation of Dolph-Chebyshev window or, for short, Dolph Window, using

!  the expression in the reference:

!    Antoniou, Andreas, 1993: Digital Filters: Analysis,

!    Design and Applications. McGraw-Hill, Inc., 689pp.

!

!  The Dolph window is optimal in the following sense:

!  For a given main-lobe width, the stop-band attenuation is minimal;

!  for a given stop-band level, the main-lobe width is minimal.

!

! Method:

!

! Modules used:    NONE

!

!------------------------------------------------------------------------------

! Parameter List:

! ---------------

INTEGER (KIND=iintegers), INTENT (IN)             ::  &

  m                   ! for dimensioning the work arrays

REAL  (KIND=ireals), INTENT (IN)                  ::  &

  deltat, taus        ! time step and cutoff period for filtering

REAL  (KIND=ireals), INTENT (OUT)                 ::  &

  window(0:2*m)       ! result

! The following variables are only used for work space

REAL  (KIND=ireals), INTENT (OUT)                 ::  &

  t(0:2*m), time(0:2*m), time2(0:2*m), w(0:2*m), w2(0:2*m)

! Local Variables:

! ----------------

INTEGER (KIND=iintegers)  :: nt, i, n, nm1, nn

REAL    (KIND=ireals)     :: zpi, zthetas, zx0, zarg, zterm1, zterm2, zrr,   &

                             zr, zdb, zsum, zsumw

!------------ End of header ---------------------------------------------------

! Begin subroutine dolph

  zpi = 4.0_ireals * ATAN(1.0_ireals)

  n = 2*m+1

  nm1 = n-1

  zthetas = 2.0_ireals*zpi*deltat/taus

  zx0 = 1.0_ireals / COS(zthetas/2.0_ireals)

  zterm1 = (zx0 + SQRT(zx0**2-1))**(REAL (N-1, ireals))

  zterm2 = (zx0 - SQRT(zx0**2-1))**(REAL (N-1, ireals))

  zrr = 0.5*(zterm1 + zterm2)

  zr = 1/zrr

  zdb = 20.0_ireals * LOG10(zr)

!------------------------------------------------------------

  DO nt = 0, M

    zsum = 1

    DO i = 1, M

      zarg = zx0 * cos(i*zpi/N)

      ! Calculate the Chebyshev polynomials

      ! Reference: Numerical Recipes, Page 184, recurrence

      !   T_n(x) = 2xT_{n-1}(x) - T_{n-2}(x) ,  n>=2.

      T(0) = 1

      T(1) = zarg

      DO nn=2,nm1

        T(nn) = 2*zarg*T(nn-1) - T(nn-2)

      ENDDO

      zterm1 = T(nm1)

      zterm2 = cos(2*nt*zpi*i/n)

      zsum   = zsum + zr*2 * zterm1 * zterm2

    ENDDO

    w(nt) = zsum / n

    TIME(nt) = nt

  ENDDO

  ! Fill in the negative-time values by symmetry.

  DO nt = 0, m

    w2(m+nt) = w(nt)

    w2(m-nt) = w(nt)

    time2(m+nt) =  time(nT)

    time2(m-nt) = -time(nT)

  ENDDO

  ! Fill up the array for return

  zsumw = 0.0_ireals

  DO nt = 0, 2*m

    zsumw = zsumw + w2(nt)

  ENDDO

  DO nt=0,2*m

    WINDOW(nt) = w2(nt)

  ENDDO

!

!

!----------------------------------------------------------

!       PRINT *, (w2(nT),    nT=0,2*M)

!----------------------------------------------------------

!

END SUBROUTINE dolph

!==============================================================================

!==============================================================================

!------------------------------------------------------------------------------

SUBROUTINE elapsed_time    (realtimedif, istat)

!------------------------------------------------------------------------------

!

! Description:

!   Returns the elapsed wall-clock time in seconds since the last call. On

!   the first call the variables are only initialized. If no system clock is

!   present, an error value of istat=1 will be returned, if the optional

!   argument istat was passed from the calling routine. 

!   realtimedif is set to 0 then.

!

! Method:

!   The intrinsic function SYSTEM_CLOCK is used, that returns the number of

!   clock counts since some system dependent event in the past (e.g. midnight

!   for a 24-hour system clock). The difference of clock counts since the last

!   call is determined and converted into seconds. The variables "lfirst"

!   and "icountsold" (see below) have to be SAVEd for the next call.

!

! Modules used:    NONE

!

!------------------------------------------------------------------------------

!

! Parameter List:

! ---------------

REAL  (KIND=ireals), INTENT (OUT)                 ::  &

      realtimedif     ! wall-clock time since the last call in seconds

                      ! (0 if no system-clock is available)

INTEGER (KIND=iintegers), INTENT (OUT), OPTIONAL  ::  &

      istat           ! optional argument for error value

! Local Variables:

! ----------------

LOGICAL, SAVE      :: lfirst = .TRUE.   ! determine whether first call or not

INTEGER, SAVE      :: icountsold        ! number of counts in the last call

INTEGER            :: icountsnew,     & ! number of counts in this call

                      ir, im            ! other arguments to SYSTEM_CLOCK

LOGICAL            :: lpres             ! if optional argument is present

!------------ End of header ---------------------------------------------------

! Begin subroutine elapsed_time

  lpres = PRESENT (istat)

  CALL SYSTEM_CLOCK ( COUNT=icountsnew, COUNT_RATE=ir, COUNT_MAX=im )

  IF ( ir /= 0 ) THEN

    ! system clock is present

    IF (lpres) THEN

      istat = 0

    ENDIF

    IF (lfirst) THEN

      ! first call: store value for the number of clock counts

      icountsold = icountsnew

      lfirst     = .FALSE.

    ELSE

      ! convert the clock counts to seconds

      IF ( icountsnew >= icountsold ) THEN

        realtimedif = ( REAL (icountsnew - icountsold, ireals) )      &

                      / REAL (ir,ireals)

      ELSE

        realtimedif = REAL (im- (icountsold-icountsnew ), ireals)     &

                      / REAL (ir, ireals)

      ENDIF

      icountsold = icountsnew

    ENDIF

  ELSE

    ! no system clock present: set error value

    realtimedif = 0.0

    IF ( lpres ) THEN

      istat = 1

    ENDIF

  ENDIF

END SUBROUTINE elapsed_time

!==============================================================================

!==============================================================================

!------------------------------------------------------------------------------

SUBROUTINE get_utc_date (ntsteps, ystartdate, dt, itype_calendar,           &

                         yactdate1, yactdate2, nactday, acthour)

!------------------------------------------------------------------------------

!

! Description:

!   This routine determines the actual date of this forecast step.

!

! Method:

!   Using the date of the forecast-start, the number of time steps 

!   already performed and the length of the time steps, the actual

!   date is calculated taking leap-years into consideration.

!   The date is given in three different formats.

!

! Modules used:    NONE

!

!------------------------------------------------------------------------------

!

! Input Parameter list:

! ---------------------

INTEGER   (KIND=iintegers), INTENT(IN)   ::                           &

  itype_calendar,   & ! for specifying the calendar used

  ntsteps             ! number of actual performed time-steps

REAL      (KIND=ireals), INTENT(IN)      ::                           &

  dt         ! time step in seconds

CHARACTER (LEN=10), INTENT(IN)           ::                           &

  ystartdate ! start date of the forecast

! Output Parameter list:

! ----------------------

CHARACTER (LEN=10), INTENT(OUT)          ::                           &

  yactdate1  ! actual date in the form   yyyymmddhh

CHARACTER (LEN=22), INTENT(OUT)          ::                           &

  yactdate2  ! actual date in the form   wd   dd.mm.yy  hh UTC

INTEGER   (KIND=iintegers), INTENT(OUT)  ::                           &

  nactday    ! day of the year

REAL      (KIND=ireals), INTENT(OUT)     ::                           &

  acthour    ! actual hour of the day

! Local variables:

INTEGER   (KIND=iintegers)   ::                                       &

  month(12), monthsum(13), ileap, iweek, iy, m,                       &

  idd, imm, iyy, ihh, iday, imonth, iyear, ihour, immhours, iyyhours, &

  iyear_hours

CHARACTER (LEN=3)            :: yweek(7)

! And for computing the amount of seconds of the whole forecast time,

! an 8-Byte INTEGER has to be used. Otherwise the computation fails after

! approx. 68 years!!

INTEGER, PARAMETER :: int_dp = KIND(1_8)

REAL(KIND=ireals)  :: zseconds

!------------ End of header ---------------------------------------------------

! Begin subroutine get_utc_date

DATA         month  / 31 ,  28 ,  31 ,  30 ,  31 ,  30 ,       &

                      31 ,  31 ,  30 ,  31 ,  30 ,  31 /

DATA         yweek  /'MON', 'TUE', 'WED', 'THU', 'FRI', 'SAT', 'SUN' /

! Statementfunction: ileap(yy) = 0:  no leap year, 

!                    ileap(yy) = 1:  leap year

! corrected version for Gregorian / Proleptic calendar

! by A. Dobler, CLM Community

  ileap (iy) = IABS( MOD(iy,  4) -   4) /   4  & ! every       4 years is a leapyear

              -IABS( MOD(iy,100) - 100) / 100  & ! every     100 years is no leapyear

              +IABS( MOD(iy,400) - 400) / 400    ! but every 400 years is a leapyear

! Divide ystartdate in day, month, year and hour

! and calculate the sums of days from the beginning of the year to the 

! end of the months

  READ ( ystartdate, '(I4,3I2)' ) iyy, imm, idd, ihh

  IF     (itype_calendar == 0) THEN

    month (2)    = 28 + ileap (iyy)

    monthsum(1) =  0

    DO m =  2 , 13

      monthsum(m) =  monthsum(m-1) + month(m-1)

    enddo

  ELSEIF (itype_calendar == 1) THEN

    monthsum(1) =  0

    DO m =  2 , 13

      monthsum(m) =  monthsum(m-1) + 30

    enddo

  ENDIF

! Determine how many hours have passed in this year

  iyyhours = (idd*24) + monthsum(imm)*24 + (ihh-24)

  iyyhours = iyyhours +                                           &

             INT (NINT (ntsteps*dt, int_dp)/3600_int_dp, iintegers)

! Take turning of the year into account

  IF     (itype_calendar == 0) THEN

    iyear_hours = 8760 + ileap(iyy)*24.0_ireals

  ELSEIF (itype_calendar == 1) THEN

    iyear_hours = 8640

  ENDIF

  IF (iyyhours < 0) THEN

    iyear    = iyy-1

    IF     (itype_calendar == 0) THEN

      iyyhours = 8760 + ileap(iyear)*24.0_ireals + iyyhours

    ELSEIF (itype_calendar == 1) THEN

      iyyhours = 8640                            + iyyhours

    ENDIF

  ELSE IF (iyyhours >= iyear_hours) THEN

    ! Take also into account if the run lasts for several years

    iyear    = iyy

    IF     (itype_calendar == 0) THEN

      iyear_hours = 8760 + ileap(iyear)*24.0_ireals

    ELSEIF (itype_calendar == 1) THEN

      iyear_hours = 8640

    ENDIF

    DO WHILE (iyyhours >= iyear_hours)

      iyyhours = iyyhours - iyear_hours

      iyear=iyear+1

      IF     (itype_calendar == 0) THEN

        iyear_hours = 8760 + ileap(iyear)*24.0_ireals

      ELSEIF (itype_calendar == 1) THEN

        iyear_hours = 8640

      ENDIF

    ENDDO

  ELSE

    iyear    =   iyy

  ENDIF

  ! calculate monthsum for actual year

  IF     (itype_calendar == 0) THEN

    month (2)    = 28 + ileap (iyear)

    monthsum(1) =  0

    DO m =  2 , 13

      monthsum(m) =  monthsum(m-1) + month(m-1)

    enddo

  ELSEIF (itype_calendar == 1) THEN

    monthsum(1) =  0

    DO m =  2 , 13

      monthsum(m) =  monthsum(m-1) + 30

    enddo

  ENDIF

! Determine the actual date from iyyhours

  m        = 1

  immhours = iyyhours

  DO WHILE (immhours >= 0)

    m        = m+1

    immhours = iyyhours - monthsum(m) * 24

  ENDDO

  imonth   = m-1

  immhours = iyyhours - monthsum(imonth)*24

  iday     = immhours/24 + 1

  ihour    = MOD(immhours,24)

!US: This was before, when 3600.0/dt was an integer value

!  acthour  = REAL (ihour, ireals) + dt/3600._ireals*  &

!                   MOD(ntsteps,INT(3600./dt+0.01)) + 0.0001_ireals

! this is the more accurate computation

  zseconds = REAL (ntsteps, ireals) * dt / 3600.0_ireals

  acthour  = REAL (ihour, ireals) +                          &

                  (zseconds - REAL(INT(zseconds, int_dp), ireals))

  ihour    = INT(acthour)

  nactday  = monthsum(imonth) + iday + INT(acthour/24. + 0.0001)

  iweek    = MOD(monthsum(imonth) + iday + (iyear-1901) + (iyear-1901)/4, 7)+1

  WRITE ( yactdate1(1:4) , '(I4.4)' ) iyear

  WRITE ( yactdate1(5:6) , '(I2.2)' ) imonth

  WRITE ( yactdate1(7:8) , '(I2.2)' ) iday

  WRITE ( yactdate1(9:10), '(I2.2)' ) ihour