Skip to content
Snippets Groups Projects
Select Git revision
  • 0ffacae75baeaef6a5d3c1a3f6788188c3bc9817
  • main default
2 results

saxon9he.jar

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    cat.f90 56.68 KiB
    ! ****************
    ! * MODUL CATMOD *
    ! ****************
    module catmod
    
    ! *********************************************
    ! *                                           *
    ! *                 CAT                       * 
    ! *                                           *
    ! *        Computer Aided Turbulence          *
    ! *                                           *
    ! *        Version  1.1    July 2016          *
    ! *                                           *
    ! *********************************************
    ! *  Theoretical Meteorology - KlimaCampus    * 
    ! *         University of Hamburg             *
    ! *********************************************
    ! *     Hartmut Borth - Edilbert Kirk         *
    ! *                                           *
    ! *           Valerio Lucarini                *
    ! *********************************************
    
    ! *********************************************
    ! * The latest version of the code and the    *
    ! * the User's Guide can be downloaded from   *
    ! * the github repository                     *
    ! *                                           *
    ! *  https://github.com/HartmutBorth/QUAD     *
    ! *********************************************
    
    
    ! *********************************************
    ! * For more details on parameters, variables *
    ! * and default values, see the User's Guide, *
    ! * which is included in the distribution of  *
    ! * CAT.                                      *
    ! *********************************************
    
    ! *****************
    ! * Model control *
    ! *****************
    character (256) :: catversion = "July 2016, Version 0.1"
    
    integer :: nshutdown = 0       ! Flag to stop program
    logical :: lrsttest  = .false. ! flag set to .true. for restart test, 
                                   ! parameter nsteps is not written 
                                   ! to restart file. lrsttest is set to
                                   ! true if a file RSTTEST is found in the
                                   ! run directory. Call <make rsttest> to
                                   ! do the restart test. To clean the run after
                                   ! a restart test call <make cleanresttest> 
    
    
    ! ***************************
    ! * Physics and mathematics *
    ! ***************************
    real(8), parameter :: pi    = 4.0d0 * atan(1.0d0)
    real(8), parameter :: twopi = pi+pi
    
    complex(8), parameter :: ci = (0.0d0,1.0d0)  ! complex unit
    
    
    ! ****************
    ! * Input/output *
    ! ****************
    
    !--- flags and switches
    logical :: lrst    = .false.   ! true if <cat_rstini> exists
    logical :: lcatnl  = .false.   ! true if <cat_namelist> exists
    
    integer :: ios_catnl  = 1      ! 0 if <cat_namelist> is readable
    
    
    !--- i/o units
    integer, parameter :: nucatnl     = 10  ! cat namelist
    integer, parameter :: nuin        = 15  ! initial conditions
    integer, parameter :: nutseri     = 20  ! time series
    integer, parameter :: nucfl       = 25  ! time series
    integer, parameter :: nugp        = 30  ! gridpoint output fields
    integer, parameter :: nusp        = 35  ! spectral output fields
    integer, parameter :: nurstini    = 40  ! restart for reading initial state
    integer, parameter :: nurstfin    = 45  ! restart for writing final state
    integer, parameter :: nudiag      = 50  ! statistics of model run
    
    integer, parameter :: nutmp       = 999 ! temporary unit
    
    
    !--- i/o file names
    character (256) :: cat_namelist  = "cat_namelist"
    character (256) :: cat_tseri     = "cat_tseri"
    character (256) :: cat_cfl       = "cat_cfl"
    character (256) :: cat_gp        = "cat_gp"
    character (256) :: cat_sp        = "cat_sp"
    character (256) :: cat_rstini    = "cat_rstini"
    character (256) :: cat_rstfin    = "cat_rstfin"
    character (256) :: cat_diag      = "cat_diag"
    
    !--- header for output of service format
    integer :: ihead(8)
    
    
    !--- codes of variables
    integer, parameter :: qcde       = 138     ! potential vorticity code
    integer, parameter :: qfrccde    = 1138    ! potential vorticity forcing code
    
    
    !--- codes of variables to be read at initialization
    integer, parameter :: ningp  = 5
    integer, parameter :: ninsp  = 5
    integer :: ingp(ningp) = &
       (/ qcde       , &
          qfrccde    , &     
           0         , &     
           0         , &     
           0           &    
       /) 
    integer :: insp(ninsp) = &
       (/  qcde      , &
           qfrccde   , &
           0         , &     
           0         , &     
           0           &
       /) 
    
    !--- codes of variables to be written to cat_gp or cat_sp 
    integer, parameter :: noutgp = 5
    integer, parameter :: noutsp = 5
    integer :: outgp(noutgp) = &
       (/ qcde   , &
           0     , &
           0     , &
           0     , &
           0       &
       /)
    integer :: outsp(noutsp) = &
       (/  0  , &
           0  , &
           0  , &
           0  , &
           0    &
       /)
    
    
    ! *******************************
    ! * Basic diagnostic parameters * 
    ! *******************************
    
    integer :: tsps          ! time steps per second (cpu-time)
    
    real    :: tmstart       ! time at start of cpu-time keeping
    real    :: tmstop        ! time at stop of cpu-time keeping
    real    :: tmrun         ! total cputime
    
    
    ! *************************************
    ! * Basic parameters of model physics *
    ! *************************************
    
    !--- non-dimensional size of fluid domain
    real(8) :: lx = twopi   ! in x-direction (scale L_x = X/2pi)
    real(8) :: ly = twopi   ! in y-direction (scale L_x = X/2pi)
    
    !--- parameters of evolution equations
    real(8) :: alpha = 0.0  ! alpha = 1/R_hat**2 [1/m**2]
                            ! with the non-dimensional Rossby
                            ! radius Ro_hat = Ro/L_x
                                
    real(8) :: beta  = 0.0  ! ambient vorticity gradient [1/m*s]
    
    !---------------------!
    ! dissipation methods !
    !---------------------!
    integer :: diss_mthd  = 1  ! dissipation method
                               ! 1: Laplacian viscosity and friction
                               !    characterized by the coefficients
                               !    sig and lam and the powers psig and
                               !    plam
                               ! 2: Laplacian viscosity and friction
                               !    characterized by the reciprocal 
                               !    damping time-scales rtsig and rtlam,
                               !    the cut-off wave-numbers ksig and klam
                               !    and the powers psig and plam 
     
    
    !----------------------------------!
    ! Laplacian viscosity and friction !
    !----------------------------------!
    
    integer, parameter :: nsig = 2 ! maximum number of terms of Laplacian 
                                       ! dissipation on small scales
     
    integer, parameter :: nlam = 2 ! maximum number of terms of Laplacian 
                                       ! dissipation on large scales
    
    !--- definition of dissipation on small scales
    real(8) :: psig (nsig) =      & ! powers of Laplacian parameterizing
           (/ 1.d0,               & ! small-scale dissipation psig > 0
              4.d0                &
           /)
    real(8) :: sig(nsig)   =      & ! coefficients of different powers
           (/ 9.765625d-05,       & ! of the Laplacian [m^2*psig/s]
              9.094947d-11        & 
           /)
    real(8) :: rtsig(nsig) =      & ! 1/time-scale of different powers
           (/ 1.d-1,              & ! of the Laplacian [1/s]
              1.d+2               & 
           /)
    integer :: ksig (nsig) =      & ! lower "cut-off" wave number of different
           (/ 320,                & ! powers of Laplacian
              100                 &
           /)
    
    !--- definition of dissipation on large scales
    real(8) :: plam (nlam) =      & ! powers of Laplacian modelling viscosity
           (/ -1.d0,              & ! on large scales plam <= 0
              -4.d0               &
           /)      
    real(8) :: lam(nlam)   =      & ! coefficients of different powers
           (/ 0.d0,               & ! of the Laplacian [m^(-2*psig)/s]
              0.d0                & 
           /)
    real(8) :: rtlam(nlam) =      & ! 1/time-scale of different powers
           (/ 0.d0,               & ! of the Laplacian [1/s]
              0.d0                & 
           /)
    integer :: klam (nlam) =      & ! lower "cut-off" wave number of different
           (/ 1,                  & ! powers of Laplacian
              1                   &
           /)
    
    !--- Random number generation
    integer, parameter  :: mxseedlen         = 8 ! maximum seed length
    integer             :: nseedlen          = 0 ! seed length
    integer             :: myseed(mxseedlen) = 0 ! seed given by namelist
    integer,allocatable :: seed(:)               ! seed defined by clock
    
    !--- Perturbation of initial conditions
    integer :: npert = 0 ! initial perturbation switch
                         ! all perturbations have zero mean-vorticity
                         !
                         ! 0 = no perturbation
                         ! 1 = white noise perturbation
                         ! 2 = white noise symmetric perturbation about the axis
                         !     y = pi (zero mean vorticity in the upper and
                         !     lower part of the fluid domain separately)
                         ! 3 = white noise anti-symmetricy perturbation about 
                         !     the axis y = pi (zero mean vorticity in the upper 
                         !     and lower power separately)
                         ! 4 = white noise symmetric perturbation about the axis
                         !     x = pi (zero mean vorticity in the upper and
                         !     lower part of the fluid domain separately)
                         ! 5 = white noise anti-symmetricy perturbation about 
                         !     the axis x = pi (zero mean vorticity in the upper 
                         !     and lower power separately)
     
    real(8) :: apert = 1.0d-5 ! amplitude of perturbation
    
    
    !--- Forcing
    real(8) :: aforc = 0.015     ! amplitude of forcing
    real(8) :: tforc = 0.001     ! memory time scale of forcing [s]
    
    integer :: in(1000,2)
    integer :: nk
    integer :: itau
    
    real(8)    :: ampcoeff
    
    complex(8) :: phi
    
    integer :: nforc = 0 ! forcing switch
                         ! 0 = no forcing
                         !
                         ! 1 = constant frocing defined by external field
                         ! 2 = spectral forcing with constant amplitudes
                         !     forcing ring and random phases
                         ! 3 = markov chain
                         ! 4 = forcing with constant wave numbers
                         !     but random uncorrelated phases (white noise)
    
    integer :: kfmin  = 0   ! min radius = sqrt(kfmin) of spectr. forcing
    integer :: kfmax  = 8   ! max radius = sqrt(kfmax) of spectr. forcing
    
    !--- Scaling factors for different parts of the evolution equation
    real(8) :: jac_scale = 1.0 ! scaling factor for jacobian
    
    
    ! ***************************************!
    ! * Basic parameters of numerical scheme !
    ! ***************************************!
    
    !--- grid in physical space
    integer :: ngx = 64    ! number of grid points (x-direction)
    integer :: ngy = 64    ! number of grid points (y-direction)
    
    integer :: nxy  = 4096                      ! ngx*ngy
    real(8) :: rnx  = 1.56250000000000000E-002  ! 1/ngx
    real(8) :: rny  = 1.56250000000000000E-002  ! 1/ngy
    real(8) :: rnxy = 2.44140625000000000E-004  ! 1/nxy
    
    
    !--- grid in spectral space
    integer :: nkx  = 21  ! ngx/3 (max. wave number in x-direction)
    integer :: nky  = 21  ! ngy/3 (max. wave number in y-direction)
    integer :: nfx  = 43  ! 2*nkx+1 (x-dimension in fourier domain)
    integer :: nfy  = 42  ! 2*nky   (y-dimension in fourier domain)
    
    
    !--- Jacobian
    integer :: jacmthd  = 1       ! approximation method of Jacobian
                                  ! 0 : no Jacobian
                                  ! 1 : divergence form (uq)_x + (vq)_y
                                  ! 2 : advection form   uq_x  + vq_y
                                  ! 3 : invariant form   
    
    !--- time integration
    integer :: nsteps    = 10000  ! number of time steps to be integrated
    integer :: tstep     = 0      ! current time step (since start of runs) 
    integer :: tstop     = 0      ! last time step of run
    
    integer :: tstp_mthd = 1      ! time stepping method
                                  ! 1 = third order Adams-Bashford
    
    integer :: nstdout   = 2500   ! time steps between messages to
                                  ! standard output (-1 no output)
    
    
    integer :: ndiag     = 500    ! time steps between test outputs into
                                  ! out_diag (-1 no test outputs) 
    
    integer :: ngp       = 100     ! time steps between output of grid point
                                   ! fields (-1 = no output)
    integer :: nsp       = 100     ! time steps between output of spectral
                                   ! fields (-1 = no output)
    integer :: ncfl      = 100     ! time steps between cfl-check 
                                   ! (-1 no cfl check)
    
    integer :: ntseri    = 100     ! time steps between time-series output
                                   ! (-1 no time-series)
     
    real(8) :: dt        = 1.d-3   ! length of time step [s]
    
    !--------------------------------------------!
    ! variables in physical/gridpoint space (GP) !
    !--------------------------------------------!
    
    !--- basic model variables
    real(8), allocatable :: gq(:,:)    ! vorticity            [1/s]
    real(8), allocatable :: gpsi(:,:)  ! stream function      [m^2/s]
    real(8), allocatable :: gu(:,:)    ! velocity x-direction [m/s]
    real(8), allocatable :: gv(:,:)    ! velocity y-direction [m/s]
    
    
    !--- jacobian (utility variables)
    
    !--- case 1
    real(8), allocatable :: guq(:,:)    ! u*q  [m/s^2]
    real(8), allocatable :: gvq(:,:)    ! v*q  [m/s^2]
    
    !--- case 2 (!! variables contain different fields depending on context)
    real(8), allocatable :: gqx (:,:) ! dq/dx [1/ms] 
    real(8), allocatable :: gqy (:,:) ! dq/dy [1/ms]
    real(8), allocatable :: gjac(:,:) ! jacobian [s^-2]
    
    !--- case 3
    real(8), allocatable :: gv2mu2(:,:) ! v^2-u^2 [m^2/s^2]
    real(8), allocatable :: guv(:,:)    ! u*v     [m^2/s^2]
    
    !--- gui
    real(4), allocatable :: ggui(:,:) ! single precision transfer
    
    
    !----------------------------------------------------------!
    ! variables in spectral space, representation as imaginary !
    ! and real part (prefix f) used for fourier transformation ! 
    ! and complex representation (prefix c) used for time      !
    ! evolution                                                !
    !----------------------------------------------------------!
    
    !--- basic model variables (format used for fourier transform)
    real(8), allocatable :: fpsi(:,:) ! stream function
    real(8), allocatable :: fu(:,:)   ! velocity x-direction
    real(8), allocatable :: fv(:,:)   ! velocity y-direction
    
    !--- jacobian (utility variables)
    
    !--- case 1
    real(8), allocatable :: fuq    (:,:) ! u*q  [m/s^2]
    real(8), allocatable :: fvq    (:,:) ! v*q  [m/s^2]
    
    !--- case 2 (!! variables contain different fields depending on context)
    real(8), allocatable :: fqx(:,:) ! u*dq/dx [s^-2] or dq/dx [1/ms]
    real(8), allocatable :: fqy(:,:) ! v*dq/dy [s^-2] or dq/dy [1/ms]
    
    !--- case 3
    real(8), allocatable :: fv2mu2 (:,:) ! v^2-u^2 [m^2/s^2]
    real(8), allocatable :: fuv    (:,:) ! u*v     [m^2/s^2]
    
    
    !--- basic model variables (format used for time propagation)
    complex(8), allocatable :: cq(:,:)     ! vorticity
    
    complex(8), allocatable :: cjac0(:,:)  ! Jacobian at time  0
    complex(8), allocatable :: cjac1(:,:)  ! Jacobian at time -1
    complex(8), allocatable :: cjac2(:,:)  ! Jacobian at time -2
    
    complex(8), allocatable :: cqfrc(:,:) ! external vorticity forcing
    
    !--- operators in Fourier Space
    integer   , allocatable :: kx(:),ky(:)
    real(8)   , allocatable :: kx2(:), ky2(:)
    real(8)   , allocatable :: k2(:,:), k2pa(:,:), rk2pa(:,:)
    real(8)   , allocatable :: kxrk2pa(:,:),kyrk2pa(:,:)
    real(8)   , allocatable :: kx2mky2(:,:), kxky(:,:)
    
    
    !--- linear time propagation (dissipation and beta-term)
    complex(8), allocatable :: cli(:,:)    ! linear time propagation
    
    
    ! *******************
    ! * External moduls *
    ! *******************
    
    !--- gui communication (guimod)
    integer :: ngui    = 1   ! global switch and parameter 
                             ! ngui > 0 GUI = on  
                             ! ngui specifies moreover the number of 
                             ! time-steps between two calls of gui_transfer 
                             ! which exchanges data between GUI and CAT
    
    integer :: nguidbg = 0   ! GUI debug mode
    
    integer :: ndatim(6) = 1 ! date/time display
    real(4) :: parc(5) = 0.0 ! timeseries display
    
    
    !--- code of predefined simulations (simmod)
    integer :: nsim  = 0 ! 0 no predefined simulation is specified 
                         ! nsim > 0 predefined simulation is run
                         ! 
                         ! available predefined simulations (experiments):
                         ! -----------------------------------------------
                         !    code       name
                         !
                         !     1    decaying_jet01
                         !    51    top-hat wind-forcing
                         !    52    top-hat wind-forcing
                         ! ----------------------------------------------- 
                         ! 
                         ! Predefined simulations are given in the dat
                         ! directory as sim_XXXX.nl with XXXX the code  
                         ! To activate a given simulation one has to copy
                         ! a given sim_XXXX.nl to sim_namelist in the run 
                         ! directory.
                         ! ----------------------------------------------- 
    
    !--- usermod (usermod)
    integer :: nuser = 0 ! 1/0 user mode is switched on/off.
                         ! In user mode it is possible to introduce 
                         ! new code into CAT (see chapter <Modifying 
                         ! CAT> in the User's Guide).
    
    
    !--- postprocessing (included in cat)
    integer :: npost = 0 ! 1/0 post processing is switched on/off
    
    
    !--- grid point fields
    real(8), allocatable :: gqxm(:)   ! mean vorticity allong x-dir [1/s]
    real(8), allocatable :: gqym(:)   ! mean vorticity allong y-dir [1/s]
    
    real(8), allocatable :: gpsixm(:) ! mean stream function allong x-dir [m^2/s]
    real(8), allocatable :: gpsiym(:) ! mean stream function allong y-dir [m^2/s]
    
    real(8), allocatable :: guxm(:)   ! mean x-velocity allong x-dir [m/s]
    real(8), allocatable :: guym(:)   ! mean x-velocity allong y-dir [m/s]
    
    real(8), allocatable :: gvxm(:)   ! mean y-velocity allong x-dir [m/s]
    real(8), allocatable :: gvym(:)   ! mean y-velocity allong y-dir [m/s]
    
    real(4), allocatable :: gguixm(:) ! single precision mean in x-dir for gui 
    real(4), allocatable :: gguiym(:) ! single precision mean in y-dir for gui 
    
    end module catmod
    
    
    ! *******
    ! * CAT *
    ! *******
    program cat
    use catmod
    implicit none
    
    call cat_prolog
    
    call cat_master
    
    call cat_epilog
    
    stop
    end program cat
    
    
    ! #####################################
    ! #      SUBROUTINES & FUNCTIONIS     #
    ! #####################################
    
    
    ! *************************
    ! * SUBROUTINE CAT_PROLOG *
    ! *************************
    subroutine cat_prolog
    use catmod
    implicit none
    
    call cpu_time(tmstart)
    inquire(file="RSTTEST",exist=lrsttest)
    call init_rst
    call init_diag
    call read_resol
    call init_pars
    call cat_alloc
    call init_ops
    if (lrst) then
       call check_rst
       call read_rst
    endif
    call cat_readnl
    call cat_open
    if (nsim > 0) call simstart
    if (nuser > 0) call userstart
    call cat_input
    call init_jacobian
    call init_lintstep
    call init_rand
    call init_forc
    call add_pert
    call init_tstep
    call init_post
    if (ngui > 0) call guistart
    
    return
    end subroutine cat_prolog
    
    
    ! *************************
    ! * SUBROUTINE READ_RESOL *
    ! *************************
    subroutine read_resol
    use catmod
    
    character (80) :: argtmp
    
    call get_command_argument(1,argtmp)
    read(argtmp,*) ngx
    
    return
    end subroutine read_resol
    
    
    ! ************************
    ! * SUBROUTINE INIT_DIAG *
    ! ************************
    subroutine init_diag
    use catmod
    implicit none
    
    open(nudiag,file=cat_diag)
    
    write(nudiag, &
     '(" *************************************************")')
    write(nudiag, &
     '(" * CAT ",a40," *")') trim(catversion)
    write(nudiag, &
       '(" *************************************************",/)')
    
    
    if (lrsttest) then
       write(nudiag, &
        '(" *************************************************")')
       write(nudiag, &
        '(" * !!! Run CAT in restart test mode !!!")')
       write(nudiag, &
        '(" *************************************************",/)')
    endif
    
    return
    end subroutine init_diag
    
    
    ! ***********************
    ! * SUBROUTINE READ_RST *
    ! ***********************
    subroutine read_rst
    use catmod
    
    write(nudiag, &
     '(/," ************************************************")')
    write(nudiag, &
     '("  Reading parameters from restart file <cat_rstini>")')
    write(nudiag, &
     '(" ************************************************",/)')
    
    
    !--- namelist parameters
    
    if (.not. lrsttest) call get_restart_iarray('nsteps',nsteps,1,1)
    call get_restart_iarray('ngp',ngp,1,1)
    call get_restart_iarray('nsp',nsp,1,1)
    call get_restart_iarray('ingp',ingp,ningp,1)
    call get_restart_iarray('insp',insp,ninsp,1)
    call get_restart_iarray('outgp',outgp,noutgp,1)
    call get_restart_iarray('outsp',outsp,noutsp,1)
    call get_restart_iarray('tstp_mthd',tstp_mthd,1,1)
    call get_restart_iarray('ncfl',ncfl,1,1)
    call get_restart_array ('dt',dt,1,1)
    call get_restart_array ('alpha',alpha,1,1)
    call get_restart_array ('beta',beta,1,1)
    call get_restart_array ('lx',lx,1,1)
    call get_restart_array ('ly',ly,1,1)
    call get_restart_array ('sig',sig,nsig,1)
    call get_restart_array ('psig',psig,nsig,1)
    call get_restart_array ('rtsig',rtsig,nsig,1)
    call get_restart_iarray('ksig',ksig,nsig,1)
    call get_restart_array ('lam',lam,nlam,1)
    call get_restart_array ('plam',plam,nlam,1)
    call get_restart_array ('rtlam',rtlam,nlam,1)
    call get_restart_iarray('klam',klam,nlam,1)
    call get_restart_iarray('diss_mthd',diss_mthd,1,1)
    call get_restart_iarray('npert',npert,1,1)
    call get_restart_array ('apert',apert,1,1)
    call get_restart_iarray('nforc',nforc,1,1)
    call get_restart_iarray('kfmin',kfmin,1,1)
    call get_restart_iarray('kfmax',kfmax,1,1)
    call get_restart_array ('aforc',aforc,1,1)
    call get_restart_array ('tforc',tforc,1,1)
    call get_restart_iarray('myseed',myseed,mxseedlen,1)
    call get_restart_iarray('ntseri',ntseri,1,1)
    call get_restart_iarray('nstdout',nstdout,1,1)
    call get_restart_iarray('jacmthd',jacmthd,1,1)
    call get_restart_iarray('ndiag',ndiag,1,1)
    call get_restart_array ('jac_scale',jac_scale,1,1)
    
    !--- additional parameters
    call get_restart_iarray('ngx',ngx,1,1)
    call get_restart_iarray('tstep',tstep,1,1)
    call get_restart_iarray('seed',seed,nseedlen,1)
    
    !--- fluid state
    call get_restart_carray('cq',cq,nkx+1,nfy+1)
    call get_restart_carray('cjac1',cjac1,nkx+1,nfy+1)
    call get_restart_carray('cjac2',cjac2,nkx+1,nfy+1)
    
    return
    end subroutine read_rst
    
    
    ! *************************
    ! * SUBROUTINE CAT_READNL *
    ! *************************
    subroutine cat_readnl
    use catmod
    implicit none
    
    
    namelist /cat_nl/ nsteps    ,ngp       ,ngui    ,ingp    ,insp       , &
                      ncfl      ,dt        ,alpha   ,beta     ,lx        , &
                      ly        ,sig       ,psig    ,rtsig    ,ksig      , &
                      lam       ,plam      ,rtlam   ,klam     ,diss_mthd , &
                      nforc     ,kfmin     ,kfmax   ,aforc    ,tforc     , &
                      myseed    ,ntseri    ,nstdout ,jacmthd  ,ndiag     , &
                      jac_scale ,nsp       ,outgp   ,outsp    ,tstp_mthd , &
                      nsim      ,nuser     ,npost   ,npert    ,apert     , &
                      nguidbg
    
    inquire(file=cat_namelist,exist=lcatnl)
    
    if (lcatnl) then
      open(nucatnl,file=cat_namelist,iostat=ios_catnl)
      read(nucatnl,cat_nl)
    endif
    
    !--- check consistency of namelist parameters   
    if (jacmthd .lt. 0 .or. jacmthd .gt. 3) then
       write(nudiag, &
       '(/," ************************************************")')
       write(nudiag,*) "jacmthd = ", jacmthd ," does not exist"
       write(nudiag,*) "use default instead (jacmthd = 1)"
       write(nudiag,'("                         ")')
       write(nudiag, &
       '(" ************************************************",/)')
       jacmthd = 1
    endif
    
    
    return
    end subroutine cat_readnl
    
    
    ! ***********************
    ! * SUBROUTINE CAT_OPEN *
    ! ***********************
    subroutine cat_open
    use catmod
    implicit none
    
    
    open(nudiag,file=cat_diag)
    
    if (ntseri .ge. 0)  open(nutseri,file=cat_tseri)
    if (ncfl .gt. 0)    open(nucfl,file=cat_cfl)
    if (ngp .gt. 0)     open(nugp,file=cat_gp,form='unformatted')
    if (nsp .gt. 0)     open(nusp,file=cat_sp,form='unformatted')
    
    return
    end subroutine cat_open
    
    
    ! ************************
    ! * SUBROUTINE INIT_PARS *
    ! ************************
    subroutine init_pars
    use catmod
    implicit none
    
    ngy  = ngx           ! restrict model to square domains
    
    nxy  = ngx * ngy     ! # of gridpoints
    nkx  = ngx / 3       ! max wavenumber in x
    nky  = ngy / 3       ! max wavenumber in y
    nfx  = nkx * 2 + 1   ! x dimension in fourier domain
    nfy  = nky * 2       ! y dimension in fourier domain
    rnx  = 1.0d0 / ngx   ! reverse of ngx
    rny  = 1.0d0 / ngy   ! reverse of ngy
    rnxy = 1.0d0 / nxy   ! reverse of gridpoint number
    
    
    write(nudiag, &
     '(" *************************************************")')
    write(nudiag, &
     '(" * Values of basic parameters                    *")')
    write(nudiag, &
     '(" *************************************************")')
    write (nudiag,*) "nkx  = ", ngx/3
    write (nudiag,*) "nky  = ", ngy/3
    write (nudiag,*) "nfx  = ", nfx
    write (nudiag,*) "nfy  = ", nfy
    write (nudiag,*) "rnx  = ", rnx
    write (nudiag,*) "rny  = ", rny
    write (nudiag,*) "rnxy = ", rnxy
    write(nudiag, &
     '(" *************************************************",/)')
    
    return
    end subroutine init_pars
    
    
    ! ************************
    ! * SUBROUTINE CAT_ALLOC *
    ! ************************
    subroutine cat_alloc
    use catmod
    implicit none
    
    integer :: j
    
    !-----------------------------
    ! 1D real and complex vectors
    !-----------------------------
    allocate(kx(0:nkx))   ; kx(:)   = [(j,j=0,nkx)]
    allocate(ky(0:nfy))   ; ky(:)   = [(j,j=0,nky),(j,j=-nky,-1)]
    allocate(kx2(0:nkx))  ; kx2(:)  = 0.0
    allocate(ky2(0:nfy))  ; ky2(:)  = 0.0
    
    
    !----------------------------
    ! 2D real and complex arrays
    !----------------------------
    
    !--- grid point space
    allocate(gq(1:ngx,1:ngy))   ; gq(:,:)   = 0.0  ! vorticity
    allocate(gu(1:ngx,1:ngy))   ; gu(:,:)   = 0.0  ! velocity in x-dir.
    allocate(gv(1:ngx,1:ngy))   ; gv(:,:)   = 0.0  ! velocity in y-dir.
    
    allocate(gpsi(1:ngx,1:ngy)) ; gpsi(:,:) = 0.0  ! stream function
    
    allocate(ggui(1:ngx,1:ngy)) ; ggui(:,:)  = 0.0  ! GUI transfer
    
    !--- spetral space 
    allocate(fu(0:nfx,0:nfy))   ; fu(:,:)   = 0.0 ! u 
    allocate(fv(0:nfx,0:nfy))   ; fv(:,:)   = 0.0 ! v 
    
    allocate(k2(0:nkx,0:nfy))      ; k2(:,:)      = 0.0 ! Laplacian
    allocate(k2pa(0:nkx,0:nfy))    ; k2pa(:,:)    = 0.0 ! modified Laplacian
    allocate(rk2pa(0:nkx,0:nfy))   ; rk2pa(:,:)   = 0.0 ! modified Laplacian^-1
    allocate(kxrk2pa(0:nkx,0:nfy)) ; kxrk2pa(:,:) = 0.0 ! q --> v
    allocate(kyrk2pa(0:nkx,0:nfy)) ; kyrk2pa(:,:) = 0.0 ! q --> u
    allocate(kx2mky2(0:nkx,0:nfy)) ; kx2mky2(:,:)  = 0.0 ! utility for Jacobian 3 
    allocate(kxky  (0:nkx,0:nfy))  ; kxky  (:,:)  = 0.0 ! utility for Jacobian 3 
    
    
    
    allocate(cli(0:nkx,0:nfy))  ; cli(:,:)   = (0.0,0.0) ! linear time propagator 
    
    allocate(cq(0:nkx,0:nfy))   ; cq(:,:)    = (0.0,0.0) ! vorticity
    allocate(cqfrc(0:nkx,0:nfy)); cqfrc(:,:) = (0.0,0.0) ! ext. vorticity force
    allocate(cjac0(0:nkx,0:nfy)); cjac0(:,:) = (0.0,0.0) ! Jacobian at time level  0
    allocate(cjac1(0:nkx,0:nfy)); cjac1(:,:) = (0.0,0.0) ! Jacobian at time level -1
    allocate(cjac2(0:nkx,0:nfy)); cjac2(:,:) = (0.0,0.0) ! Jacobian at time level -2
    
    
    !allocate(aptmp(0:nkx,0:nfy)); aptmp(:,:) = 0.0 ! amplitude of complex field
    !allocate(patmp(0:nkx,0:nfy)); patmp(:,:) = 0.0 ! phase of complex field
    
    return
    end subroutine cat_alloc
    
    
    ! ****************************
    ! * SUBROUTINE INIT_JACOBIAN *
    ! ****************************
    subroutine init_jacobian
    use catmod
    implicit none
    select case(jacmthd)
    case(0)
       print *, "init_jacobian case 0"
    case(1)
       !--- allocate fields used by jacobian
    
       !--- grid point fields
       allocate(guq(1:ngx,1:ngy)) ; guq(:,:) = 0.0  ! u*q
       allocate(gvq(1:ngx,1:ngy)) ; gvq(:,:) = 0.0  ! v*q
    
       !--- spectral fields  
       allocate(fuq(0:nfx,0:nfy)) ; fuq(:,:) = 0.0  ! u*q
       allocate(fvq(0:nfx,0:nfy)) ; fvq(:,:) = 0.0  ! v*q
    case(2)
       !--- allocate fields used by jacobian
    
       !--- grid point fields
       allocate(gqx(1:ngx,1:ngy))  ; gqx(:,:)  = 0.0  ! qx
       allocate(gqy(1:ngx,1:ngy))  ; gqy(:,:)  = 0.0  ! qy
       allocate(gjac(1:ngx,1:ngy)) ; gjac(:,:) = 0.0  ! jacobian in grid-space
    
       !--- spectral fields  
       allocate(fqx(0:nfx,0:nfy)) ; fqx(:,:) = 0.0  ! u*qx or qx
       allocate(fqy(0:nfx,0:nfy)) ; fqy(:,:) = 0.0  ! v*qy or qy
    case(3)
       !--- allocate fields used by jacobian
       !--- grid point fields 
       allocate(gv2mu2(1:ngx,1:ngy)) ; gv2mu2(:,:) = 0.0  ! v^2-u^2 
       allocate(guv   (1:ngx,1:ngy)) ; guv   (:,:) = 0.0  ! u*v 
       !--- spetral fields  
       allocate(fv2mu2(0:nfx,0:nfy)) ; fv2mu2(:,:) = 0.0  ! v^2-u^2
       allocate(fuv   (0:nfx,0:nfy)) ; fuv   (:,:) = 0.0  ! u*v
       print *, "init_jacobian case 3"
    end select
    
    return
    end subroutine init_jacobian
    
    
    ! ****************************
    ! * SUBROUTINE STOP_JACOBIAN *
    ! ****************************
    subroutine stop_jacobian
    use catmod
    implicit none
    
    select case (jacmthd)
    
    case (1)
       !--- deallocate fields used by jacobian
    
       !--- grid point fields 
       deallocate(guq)
       deallocate(gvq)
    
       !--- spetral fields 
       deallocate(fuq)
       deallocate(fvq)
    
    case (2)
       !--- deallocate fields used by jacobian
    
       !--- grid point fields 
       deallocate(gqx)
       deallocate(gqy)
       deallocate(gjac)
    
       !--- spetral fields 
       deallocate(fqx)
       deallocate(fqy)
    
    case (3)
       !--- deallocate fields used by jacobian
    
       !--- grid point fields 
       deallocate(gv2mu2)
       deallocate(guv)
    
       !--- spetral fields 
       deallocate(fv2mu2)
       deallocate(fuv)
    
    end select
    
    return
    end subroutine stop_jacobian
    
    
    ! ************************
    ! * SUBROUTINE CAT_INPUT *
    ! ************************
    subroutine cat_input
    use catmod
    implicit none
    
    logical         :: lexist
    integer         :: kcode,jj,kk,kkmax
    real(8)         :: gtmp(1:ngx,1:ngy)
    real(8)         :: ftmp(0:nfx,0:nfy) ! temporary field for i/o
    
    character(2)    :: gtp
    character(256)  :: fname
    
    integer, parameter      :: ngtp         = 2
    character(2), parameter :: gtp_ar(ngtp) = (/ "GP" , "SP" /)
    
    
    if (.not. (tstep .eq. 0) ) return
    
    do jj = 1,ngtp
       gtp = gtp_ar(jj)
       select case (gtp)
       case ("GP")
          kkmax = ningp
       case ("SP")
          kkmax = ninsp
       end select
       do kk = 1,kkmax
          select case (gtp)
          case ("GP")
             kcode = ingp(kk)
          case ("SP")
             kcode = insp(kk)
          end select
          if (kcode .gt. 0) then
             call checkvar(kcode,gtp,lexist)
             if (lexist) then
                select case (kcode)
                case (qcde)
                   select case (gtp)
                   case ("GP")
                      call read_gp(kcode,gtp,gq)
                      call grid_to_fourier(gq,cq,nfx,nfy,ngx,ngy)
                   case ("SP")
                      call read_sp(kcode,gtp,ftmp)
                      call f2c(ftmp,cq)
                   end select
                   cq(0,0) = (0.0,0.0)
                case (qfrccde)
                   select case (gtp)
                   case ("GP")
                      call read_gp(kcode,gtp,gtmp)
                      call grid_to_fourier(gtmp,cqfrc,nfx,nfy,ngx,ngy)
                   case ("SP")
                      call read_sp(kcode,gtp,ftmp)
                      call f2c(ftmp,cqfrc)
                   end select
                   cqfrc(0,0) = (0.0,0.0)
                end select
             else
                call cat_fname(gtp,kcode,fname)
                write(nudiag, &
                '(" *************************************************")')
                write(nudiag, &
                '(" *   File ", a ," not found, use default")') trim(fname)
                write(nudiag, &
                '(" *************************************************",/)')
             endif
          endif
       enddo
    enddo
    
    return
    end subroutine cat_input
    
    
    ! ************************
    ! * SUBROUTINE CAT_FNAME *
    ! ************************
    subroutine cat_fname(gtp,kcode,fname)
    use catmod
    implicit none
    
    character(2)   :: gtp
    character(256) :: fname
    integer :: kcode
    
    if (ngx < 100) then
      write(fname,'(a2,i2.2,"_var",i4.4,".srv")') gtp,ngx,kcode
    elseif (ngx < 1000) then
      write(fname,'(a2,i3.3,"_var",i4.4,".srv")') gtp,ngx,kcode
    else
      write(fname,'(a2,i4.4,"_var",i4.4,".srv")') gtp,ngx,kcode
    endif
    
    fname = trim(fname)
    
    return
    end subroutine cat_fname
    
    
    ! ***********************
    ! * SUBROUTINE CHECKVAR *
    ! ***********************
    subroutine checkvar(kcode,gtp,lexist)
    use catmod
    
    logical :: lexist
    character(2)   :: gtp
    character(256) :: fname
    integer :: kcode
    
    call cat_fname(gtp,kcode,fname)
    inquire(file=trim(fname),exist=lexist)
    
    return
    end subroutine checkvar
    
    
    ! **********************
    ! * SUBROUTINE READ_GP *
    ! **********************
    subroutine read_gp(kcode,gtp,gpvar)
    use catmod
    implicit none
    
    character(2)   :: gtp
    character(256) :: fname
    integer        :: kcode
    real(8)        :: gpvar(1:ngx,1:ngy)
    
    call cat_fname(gtp,kcode,fname)
    
    open(nuin,file=fname,form='unformatted')
    
    write(nudiag, &
    '(" *************************************************")')
    write(nudiag,'(" * Reading var",i4.4 " from file " a)') &
          kcode,trim(fname)
    write(nudiag, &
    '(" *************************************************",/)')
    
    read (nuin) ihead
    read (nuin) gpvar(:,:)
    close(nuin)
    
    return
    end subroutine read_gp
    
    
    ! **********************
    ! * SUBROUTINE READ_SP *
    ! **********************
    subroutine read_sp(kcode,gtp,spvar)
    use catmod
    implicit none
    
    character(2)   :: gtp
    character(256) :: fname
    integer        :: kcode
    real(8)        :: spvar(0:nfx,0:nfy)
    
    call cat_fname(gtp,kcode,fname)
    
    open(nuin,file=fname,form='unformatted')
    
    write(nudiag, &
    '(" *************************************************")')
    write(nudiag,'(" * Reading var",i4.4 " from file " a)') &
          kcode,trim(fname)
    write(nudiag, &
    '(" *************************************************",/)')
    
    read (nuin) ihead
    read (nuin) spvar(:,:)
    close(nuin)
    
    return
    end subroutine read_sp
    
    
    ! ******************
    ! * SUBROUTINE F2C *
    ! ******************
    subroutine f2c(fvar,cvar)
    use catmod
    implicit none
    
    integer     :: i,j
    real(8)     :: fvar(0:nfx,0:nfy)
    complex(8)  :: cvar(0:nkx,0:nfy)
    
    do j = 0, nfy
       do i = 0, nkx
          cvar(i,j) = cmplx(fvar(i+i,j),fvar(i+i+1,j))
       enddo
    enddo
    
    return
    end subroutine f2c
    
    
    ! ******************
    ! * SUBROUTINE C2F *
    ! ******************
    subroutine c2f(cvar,fvar)
    use catmod
    implicit none
    
    integer     :: i,j
    real(8)     :: fvar(0:nfx,0:nfy)
    complex(8)  :: cvar(0:nkx,0:nfy)
    
    
    
    do j = 0, nfy
       do i = 0, nkx
          fvar(i+i  ,j) = real (cvar(i,j))
          fvar(i+i+1,j) = aimag(cvar(i,j))
       enddo
    enddo
    
    return
    end subroutine c2f
    
    
    ! ****************************
    ! * SUBROUTINE INIT_LINTSTEP *
    ! ****************************
    subroutine init_lintstep
    use catmod
    implicit none
    
    real(8)    :: diss_sig,diss_lam,b_term
    complex(8) :: arg
    integer    :: jx,jy,n
    
    !--------------------------------------------------------------
    ! Time propagator of linear part of differential equation
    !
    ! cli = exp(arg)                      with
    !
    ! arg = -dt*[diss_sig + diss_lam + beta_term]
    !
    !  with 
    !     diss_sig  = -sum_j=1,nsig sig(j)*k^(2*psig(j))
    !     diss_lam  = -sum_j=1,nlam lam(j)*k^(2*plam(j))
    !     b_term    = i*beta*kx/(k^2 + alpha)
    !--------------------------------------------------------------
    
    do jy = 0, nfy
       do jx = 0, nkx
          diss_sig = 0.0
          diss_lam = 0.0
          do n = 1,nsig
             diss_sig = diss_sig - sig(n)*k2(jx,jy)**psig(n)
          enddo
          do n = 1,nlam
             diss_lam = diss_lam - lam(n)*k2(jx,jy)**plam(n)
          enddo
          b_term    = beta * kxrk2pa(jx,jy)
          arg       = dt * cmplx(diss_sig+diss_lam,b_term)
          cli(jx,jy)  = exp(arg)
       enddo
    enddo
    
    cli(0,0) = (0.0,0.0)
    
    return
    end subroutine init_lintstep
    
    
    ! ************************
    ! * SUBROUTINE INIT_RAND *
    ! ************************
    subroutine init_rand
    use catmod
    implicit none
    
    integer :: k, clock
    
    call random_seed(size=nseedlen)
    allocate(seed(nseedlen))
    
    if (myseed(1) /= 0) then
       seed(:) = 0
       k = nseedlen
       if (k .gt. mxseedlen) k = mxseedlen
       seed(1:k) = myseed(1:k)
    else
       call system_clock(count=clock)
       seed(:) = clock + 37 * (/(k,k=1,nseedlen)/)
    endif
    
    call random_seed(put=seed)
    
    return
    end subroutine init_rand
    
    
    ! ************************
    ! * SUBROUTINE INIT_FORC *
    ! ************************
    subroutine init_forc
    use catmod
    implicit none
    
    real(8) :: k2tmp,ent
    integer :: jx,jy
    
    if (nforc .eq. 0) return
    
    !---     
    ent=0.d0
    do jy = 0, nfy
      do jx = 0, nkx
    !   k2 = kx2(i) + ky2(j)
        k2tmp = k2(jx,jy)
        if (jy.ne.0.or.jx.ne.0) then
          if (sqrt(k2tmp).ge.kfmin.and.sqrt(k2tmp).le.kfmax) then
            nk = nk + 1
            in(nk,1) = jx 
            in(nk,2) = jy
    !       ent=ent+k2+alpha
            ent=ent+k2pa(jy,jx)
          endif
        endif
      enddo
    enddo
    
    ampcoeff=sqrt(aforc*dt/ent)*nxy
    
    if (tforc.le.dt) then
      itau=1
      else
        itau=tforc/dt
    endif
    
    print *, "the forcing ring contains ",nk," wavevectors."
    return
    
    
    end subroutine init_forc
    
    
    ! *************************
    ! * SUBROUTINE INIT_FORC1 *
    ! *************************
    subroutine init_forc1
    use catmod
    implicit none
    
    if (nforc .eq. 0) return
    
    
    
    !--- introduce default forcing in init_forc1
    !--- can be overwritten by input file
    
    !--- introduce markovian time-scale
    
    !--- forcing is specified by input files
    
    !--- examples are produced in simmod.f90 
    
    !--- produce masks
    
    !--- normalize forcing
    
    
    return
    end subroutine init_forc1
    
    
    ! ***********************
    ! * SUBROUTINE ADD_PERT *
    ! ***********************
    subroutine add_pert
    use catmod
    implicit none
    
    
    integer :: jx,jy                ! loop index
    real(8) :: anglex,delx          ! angle and grid distance in x-direction
    real(8) :: angley,dely          ! angle and grid distance in y-direction
    real(8) :: gqpert(1:ngx,1:ngy)  ! vorticity perturbation [1/s]
    real(8) :: xtmp(1:ngx)          ! temporary vector (x-direction)
    real(8) :: ytmp(1:ngy)          ! temporary vector (y-direction)
    
    complex(8) :: cqpert(0:nkx,0:nfy)  ! spectral vorticiy perturbation [1/s]
    
    if ( .not. (tstep .eq. 0) .and. npert .eq. 0) return
    
    delx = 4.0*asin(1.0)/ngx
    dely = 4.0*asin(1.0)/ngy
    
    select case(npert)
       !--- zero-mean grid point white noise perturbations
       case(1)
          !--- without symmetry
          call random_number(gqpert)
          gqpert = 2.0*apert*(gqpert - sum(gqpert)*rnxy)
          call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
          cqpert(0,0) = (0.0d0,0.0d0)
          cq = cq + cqpert
       case(2)
          !--- symmetric about channel center
          do jy = 1,ngy/2
             call random_number(xtmp)
             xtmp = 2.0*apert*(xtmp - sum(xtmp)*rnx)
             gqpert(:,jy)       = xtmp(:)   
             gqpert(:,ngy+1-jy) = xtmp(:)
          enddo
          call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
          cqpert(0,0) = (0.0,0.0)
          cq = cq + cqpert
       case(3)
          !--- anti-symmetric about channel center
          do jy = 1,ngy/2
             call random_number(xtmp)
             xtmp = 2.0*apert*(xtmp - sum(xtmp)*rnx)
             gqpert(:,jy)       =  xtmp(:)   
             gqpert(:,ngy+1-jy) = -xtmp(:)
          enddo
          call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
          cqpert(0,0) = (0.0,0.0)
          cq = cq + cqpert
       case(4)
          !--- symmetric about x = pi 
          do jx = 1,ngx/2
             call random_number(ytmp)
             ytmp = 2.0*apert*(ytmp - sum(ytmp)*rny)
             gqpert(jx,:)       = ytmp(:)   
             gqpert(ngx+1-jx,:) = ytmp(:)
          enddo
          call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
          cqpert(0,0) = (0.0,0.0)
          cq = cq + cqpert
       case(5)
          !--- anti-symmetric about x = pi 
    !     do jx = 1,ngx/2
    !        call random_number(ytmp)
    !        ytmp = 2.0*apert*(ytmp - sum(ytmp)*rny)
    !        gqpert(jx,:)       =  ytmp(:)   
    !        gqpert(ngx+1-jx,:) = -ytmp(:)
    !     enddo
          do jy = 1,ngy/2
             call random_number(xtmp)
             xtmp = 2.0*apert*(xtmp - sum(xtmp)*rnx)
             gqpert(:,jy)       =  xtmp(:)   
             gqpert(:,ngy+1-jy) = -xtmp(:)
          enddo
          gqpert = transpose(gqpert)
          call grid_to_fourier(gqpert,cqpert,nfx,nfy,ngx,ngy)
          cqpert(0,0) = (0.0,0.0)
          cq = cq + cqpert
       case(6)
          cqpert = (0.0d0,0.0d0)
          cqpert(1,nky/2) = 0.5d0*apert*cmplx(1.0d0,0.0d0)
          cqpert(nkx/2+1,1) = 0.5d0*apert*cmplx(1.0d0,0.0d0)
          cq = cq + cqpert
    end select
    
    return
    end subroutine add_pert
    
    
    ! *********************
    ! * SUBROUTINE CQ2GUV *
    ! *********************
    subroutine cq2guv
    use catmod
    implicit none
    
    call cq2fuv
    call fourier_to_grid(fu,gu,nfx,nfy,ngx,ngy)
    call fourier_to_grid(fv,gv,nfx,nfy,ngx,ngy)
    
    return
    end subroutine cq2guv
    
    
    ! **********************
    ! * SUBROUTINE CQ2GQUV *
    ! **********************
    subroutine cq2gquv
    use catmod
    implicit none
    
    call cq2fuv
    call fourier_to_grid(cq,gq,nfx,nfy,ngx,ngy)
    call fourier_to_grid(fu,gu,nfx,nfy,ngx,ngy)
    call fourier_to_grid(fv,gv,nfx,nfy,ngx,ngy)
    
    return
    end subroutine cq2gquv
    
    ! **********************
    ! * SUBROUTINE CQ2GQUV *
    ! **********************
    subroutine cq2gqxqy
    use catmod
    implicit none
    
    call cq2fqxqy
    call fourier_to_grid(fqx,gqx,nfx,nfy,ngx,ngy)
    call fourier_to_grid(fqy,gqy,nfx,nfy,ngx,ngy)
    
    return
    end subroutine cq2gqxqy
    
    ! *************************
    ! * SUBROUTINE CAT_WRTOUT *
    ! *************************
    subroutine cat_wrtout
    use catmod
    implicit none
    
    real(8)  :: ftmp(0:nfx,0:nfy) ! temporary field for i/o
    
    
    integer :: kk,kcode
    
    if((ngp .gt. 0) .and. mod(tstep,ngp).eq.0)  then
       do kk = 1, noutgp
          kcode = outgp(kk)
          select case (kcode)
          case (qcde)
             call cat_wrtgp(nugp,gq,qcde,0)
          end select
       enddo 
    endif
    
    if((nsp .gt. 0) .and. mod(tstep,nsp).eq.0)  then
       do kk = 1, noutsp
          kcode = outsp(kk)
          select case (kcode)
          case (qcde)
             call c2f(cq,ftmp)
             call cat_wrtsp(nusp,ftmp,qcde,0)
          end select
       enddo
    endif
    
    if((ntseri .gt. 0) .and. mod(tstep,ntseri).eq.0) call write_tseri
    
    if((ncfl .gt. 0) .and. mod(tstep,ncfl).eq.0)   call write_cfl
    
    return
    end subroutine cat_wrtout
    
    
    ! ***************************
    ! * SUBROUTINE GUI_TRANSFER *
    ! ***************************
    subroutine gui_transfer
    use catmod
    implicit none
    
    ggui(:,:) = gq(:,:) ! double precision -> single
    call guiput("GQ" // char(0), ggui, ngx, ngy, 1)
    
    if (npost > 0) then
       !--- zonal means
       gguixm(:) = gqxm(:)   ! double -> single
       call guiput("GQXM" // char(0), gguixm, ngy, 1, 1)   ! q 
    
       gguixm(:) = gpsixm(:) ! double -> single
       call guiput("GPSIXM" // char(0), gguixm, ngy, 1, 1) ! psi 
    
       gguixm(:) = guxm(:) ! double -> single
       call guiput("GUXM" // char(0), gguixm, ngy, 1, 1)   ! u
    
       gguixm(:) = gvxm(:) ! double -> single
       call guiput("GVXM" // char(0), gguixm, ngy, 1, 1)   ! v
    
    
       !--- meridional means
       gguiym(:) = gqym(:)   ! double -> single
       call guiput("GQYM" // char(0), gguiym, ngx, 1, 1) ! q 
    
       gguiym(:) = gpsiym(:) ! double -> single
       call guiput("GPSIYM" // char(0), gguiym, ngx, 1, 1) ! psi 
    
       gguiym(:) = guym(:) ! double -> single
       call guiput("GUYM" // char(0), gguiym, ngx, 1, 1)   ! u
    
       gguiym(:) = gvym(:) ! double -> single
       call guiput("GVYM" // char(0), gguiym, 1, ngx, 1)   ! v
    
    endif
    
    return
    end subroutine gui_transfer
    
    
    ! ************************
    ! * SUBROUTINE CAT_WRTGP *
    ! ************************
    subroutine cat_wrtgp(ku,gpvar,kcode,klev)
    use catmod
    implicit none
    
    integer :: ku,kcode,klev
    integer :: yy,mo,dd,hh,mm,ii
    real(8) :: gpvar(ngx,ngy)
    
    ! Build a header for service format
    
    mm = mod(tstep,60)
    ii = tstep / 60
    hh = mod(ii,24)
    ii = ii / 24
    dd = mod(ii,30) + 1
    ii = ii / 30
    mo = mod(ii,12) + 1
    yy = ii / 12
    
    ihead(1) = kcode
    ihead(2) = klev
    ihead(3) = dd + 100 * mo + 10000 * yy
    ihead(4) = mm + 100 * hh
    ihead(5) = ngx
    ihead(6) = ngy
    ihead(7) = 0
    ihead(8) = 0
    
    write (ku) ihead
    write (ku) gpvar(:,:)
    
    return
    end subroutine cat_wrtgp
    
    
    ! ************************
    ! * SUBROUTINE CAT_WRTSP *
    ! ************************
    subroutine cat_wrtsp(ku,spvar,kcode,klev)
    use catmod
    implicit none
    
    integer :: ku,kcode,klev
    integer :: yy,mo,dd,hh,mm,ii
    real(8) :: spvar(0:nfx,0:nfy)
    
    mm = mod(tstep,60)
    ii = tstep / 60
    hh = mod(ii,24)
    ii = ii / 24
    dd = mod(ii,30) + 1
    ii = ii / 30
    mo = mod(ii,12) + 1
    yy = ii / 12
    
    ihead(1) = kcode
    ihead(2) = klev
    ihead(3) = dd + 100 * mo + 10000 * yy
    ihead(4) = mm + 100 * hh
    ihead(5) = nfx+1
    ihead(6) = nfy+1
    ihead(7) = 0
    ihead(8) = 0
    
    write (ku) ihead
    write (ku) spvar(:,:)
    
    return
    end subroutine cat_wrtsp
    
    
    ! *************************
    ! * SUBROUTINE INIT_TSTEP *
    ! *************************
    subroutine init_tstep
    use catmod
    implicit none
    
    real(8) :: dt2 
    
    
    dt2 = dt/2
    
    !--- determine tstop and return if tstep > 0
    if (tstep .eq. 0) then
       tstop = tstep + nsteps      
    else
       tstop = tstep - 1 + nsteps   
       return
    endif
    
    select case (tstp_mthd)
    case (1)
       call cq2gquv
       call jacobian
       cjac1(:,:) = cjac0(:,:)
       cjac2(:,:) = cjac1(:,:)
    
       call cat_wrtout
       !--- euler-step with dt/2
       cq(:,:) = cli(:,:)*(cq(:,:)+dt2*cjac0(:,:))
    
       call cq2gquv
       call jacobian
       !--- euler-step with dt
       cq(:,:) = cli(:,:)*(cq(:,:)+dt*cjac0(:,:))
       tstep=tstep+1
    
       call cq2gquv
       call jacobian
       call cat_wrtout
    
       !--- adams-bashford method 2nd order
       cq(:,:) = cli(:,:) * (cq(:,:) + dt2 * (3.0 * cjac0(:,:) -            &
                 cli(:,:) * cjac1(:,:)))
       cq(0,0) = (0.0,0.0)
    
       if (nforc .ge. 1) call add_forc
    
       cjac1(:,:) = cjac0(:,:)
    
       tstep=tstep+1
    end select
    
    return
    end subroutine init_tstep
    
    ! ************************
    ! * SUBROUTINE INIT_POST *
    ! ************************
    subroutine init_post
    use catmod
    implicit none
    
    allocate(gqxm(1:ngy))   ; gqxm(:)   = 0.0 ! q-mean in x-dir [1/s]
    allocate(gqym(1:ngx))   ; gqym(:)   = 0.0 ! q-mean in y-dir [1/s]
    
    allocate(gpsixm(1:ngy)) ; gpsixm(:) = 0.0 ! psi-mean in x-dir [m^2/s]
    allocate(gpsiym(1:ngx)) ; gpsiym(:) = 0.0 ! psi-mean in y-dir [m^2/s]
    
    allocate(guxm(1:ngy))   ; guxm(:)   = 0.0 ! u-mean in x-dir [m/s]
    allocate(guym(1:ngx))   ; guym(:)   = 0.0 ! u-mean in y-dir [m/s]
    
    allocate(gvxm(1:ngy))   ; gvxm(:)   = 0.0 ! v-mean in x-dir [m/s]
    allocate(gvym(1:ngx))   ; gvym(:)   = 0.0 ! v-mean in y-dir [m/s]
    
    allocate(gguixm(1:ngy)) ; gguixm(:) = 0.0 ! single precision mean in x-dir  
    allocate(gguiym(1:ngx)) ; gguiym(:) = 0.0 ! single precision mean in y-dir  
    
    return
    end subroutine init_post
    
    
    ! ************************
    ! * SUBROUTINE STOP_POST *
    ! ************************
    subroutine stop_post
    use catmod
    implicit none
    
    if (npost > 0) then
       deallocate(gqxm) 
       deallocate(gqym) 
    
       deallocate(gpsixm)
       deallocate(gpsiym) 
    
       deallocate(guxm)
       deallocate(guym)
    
       deallocate(gvxm)
       deallocate(gvym)
    
       deallocate(gguixm)
       deallocate(gguiym)
    endif
    
    return
    end subroutine stop_post
    
    
    ! *************************
    ! * SUBROUTINE CAT_MASTER *
    ! *************************
    subroutine cat_master
    use catmod
    implicit none
    
    do while (tstep <= tstop)
       call cq2gquv
       call jacobian
       call post
       call cat_wrtout
       if (ngui > 0 .and. mod(tstep,ngui) == 0) then
          call gui_transfer
          call guistep_cat
       endif
       call step
       if (nforc .ge. 1) call add_forc
       if (nsim > 0) call simstep
       if (nuser > 0) call userstep
       tstep = tstep + 1
       if (nstdout.gt.0 .and. ngui == 0 .and. mod(tstep,nstdout) == 0) then
          write(*,*)' time step ',tstep
       endif
       if (nshutdown > 0) return
    enddo
    
    return
    end subroutine cat_master
    
    
    ! *******************
    ! * SUBROUTINE POST *
    ! *******************
    subroutine post
    use catmod
    implicit none
    
    !--- calculate zonal means
    gqxm   = sum(gq  ,1)*rnx
    gpsixm = sum(gpsi,1)*rnx
    guxm   = sum(gq  ,1)*rnx
    gvxm   = sum(gq  ,1)*rnx
    
    !--- calculate meridional means
    gqym   = sum(gq  ,2)*rny
    gpsiym = sum(gpsi,2)*rny
    guym   = sum(gq  ,2)*rny
    gvym   = sum(gq  ,2)*rny
    
    return
    end subroutine post
    
    
    ! *******************
    ! * SUBROUTINE STEP *
    ! *******************
    subroutine step
    use catmod
    implicit none
    
    real(8) :: c,c1
    
    
    select case (tstp_mthd)
    case (1)
       c  = dt/12.0
       c1 = 23.0 * c
    
       !--- adams-bashford 3rd order
       cq(:,:) = cli(:,:) * (cq(:,:) + c1 * cjac0(:,:) + c * cli(:,:) *  &
                 (-16.0 * cjac1(:,:) + 5.0 * cli(:,:)*cjac2(:,:)))
       cq(0,0) = (0.0,0.0)
    
       !--- shift time-levels (pointers are faster)
       cjac2(:,:) = cjac1(:,:)
       cjac1(:,:) = cjac0(:,:)
    
    end select
    
    return
    end subroutine step
    
    
    ! *************************
    ! * SUBROUTINE CAT_EPILOG *
    ! *************************
    subroutine cat_epilog
    use catmod
    implicit none
    
    if (ngui > 0) call guistop
    
    call write_rst
    
    call cpu_time(tmstop)
    tmrun = tmstop - tmstart
    tsps    = nint(nsteps / tmrun)
    
    write(nudiag, &
     '(" *************************************************")')
    write(nudiag,'("  Total time in seconds: ",f15.2)') tmrun
    write(nudiag,'("  Time steps per second: ",i12)') tsps
    write(nudiag, &
     '(" *************************************************")')
    
    call stop_jacobian
    call stop_post
    if (nsim > 0) call simstop
    if (nuser > 0) call userstop
    call close_files
    
    return
    end subroutine cat_epilog
    
    
    ! ************************
    ! * SUBROUTINE WRITE_RST *
    ! ************************
    subroutine write_rst
    use catmod
    implicit none
    
    
    !--- namelist parameters
    if (.not. lrsttest) call put_restart_iarray('nsteps',nsteps,1,1)
    call put_restart_iarray('ngp',ngp,1,1)
    call put_restart_iarray('nsp',nsp,1,1)
    call put_restart_iarray('ingp',ingp,ningp,1)
    call put_restart_iarray('insp',insp,ninsp,1)
    call put_restart_iarray('outgp',outgp,noutgp,1)
    call put_restart_iarray('outsp',outsp,noutsp,1)
    call put_restart_iarray('tstp_mthd',tstp_mthd,1,1)
    call put_restart_iarray('ncfl',ncfl,1,1)
    call put_restart_array ('dt',dt,1,1)
    call put_restart_array ('alpha',alpha,1,1)
    call put_restart_array ('beta',beta,1,1)
    call put_restart_array ('lx',lx,1,1)
    call put_restart_array ('ly',ly,1,1)
    call put_restart_array ('sig',sig,nsig,1)
    call put_restart_array ('psig',psig,nsig,1)
    call put_restart_array ('rtsig',rtsig,nsig,1)
    call put_restart_iarray('ksig',ksig,nsig,1)
    call put_restart_array ('lam',lam,nlam,1)
    call put_restart_array ('plam',plam,nlam,1)
    call put_restart_array ('rtlam',rtlam,nlam,1)
    call put_restart_iarray('klam',klam,nlam,1)
    call put_restart_iarray('diss_mthd',diss_mthd,1,1)
    call put_restart_iarray('npert',npert,1,1)
    call put_restart_array ('apert',apert,1,1)
    call put_restart_iarray('nforc',nforc,1,1)
    call put_restart_iarray('kfmin',kfmin,1,1)
    call put_restart_iarray('kfmax',kfmax,1,1) 
    call put_restart_array ('aforc',aforc,1,1)
    call put_restart_array ('tforc',tforc,1,1)
    call put_restart_iarray('myseed',myseed,mxseedlen,1)
    call put_restart_iarray('ntseri',ntseri,1,1)
    call put_restart_iarray('nstdout',nstdout,1,1)
    call put_restart_iarray('jacmthd',jacmthd,1,1)
    call put_restart_iarray('ndiag',ndiag,1,1)
    call put_restart_array ('jac_scale',jac_scale,1,1)
    
    
    !--- additional parameters
    call put_restart_iarray('ngx',ngx,1,1)
    call put_restart_iarray('tstep',tstep,1,1)
    
    call random_seed(get=seed)
    call put_restart_iarray('seed',seed,nseedlen,1)
    
    !--- fluid state
    call put_restart_carray('cq',cq,nkx+1,nfy+1)
    call put_restart_carray('cjac1',cjac1,nkx+1,nfy+1)
    call put_restart_carray('cjac2',cjac2,nkx+1,nfy+1)
    
    return
    end subroutine write_rst
    
    
    ! **************************
    ! * SUBROUTINE CLOSE_FILES *
    ! **************************
    subroutine close_files
    use catmod
    
    if (lrst)            close(nurstini)
    if (lcatnl)          close(nucatnl) 
    close(nurstfin)
    close(nudiag)
    if (ntseri .gt. 0)   close(nutseri)
    if (ncfl .gt. 0)     close(nucfl)
    if (ngp .gt. 0)      close(nugp)
    if (nsp .gt. 0)      close(nusp)
    
    return
    end subroutine close_files
    
    
    ! ***********************
    ! * SUBROUTINE JACOBIAN *
    ! ***********************
    subroutine jacobian
    use catmod
    implicit none
    
    integer    :: jx
    
    select case (jacmthd)
    case (0)
       !--- Jacobian switched off
       cjac0(:,:) = cmplx(0.0d0,0.0d0)
    
    case (1)
       !--- Jacobian in divergence form (second representation in CAT-UG)
       guq = gu*gq
       gvq = gv*gq
    
       call grid_to_fourier(guq,fuq,nfx,nfy,ngx,ngy)
       call grid_to_fourier(gvq,fvq,nfx,nfy,ngx,ngy)
    
       !--- (!!! Jacobian has different sign than in CAT-UG !!!)
       do jx = 0, nkx
          cjac0(jx,:) = & 
           cmplx( kx(jx)*fuq(jx+jx+1,:)+ky(:)*fvq(jx+jx+1,:), &
                 -kx(jx)*fuq(jx+jx  ,:)-ky(:)*fvq(jx+jx  ,:))
       enddo
    
    case (2)
       !--- Jacobian in advection form (first representation in CAT-UG)
       call cq2gqxqy
       gjac = gu*gqx + gv*gqy
       call grid_to_fourier(gjac,cjac0,nfx,nfy,ngx,ngy)
    
    case (3)
       !--- Jacobian in invariant form (third representation in CAT-UG)
       gv2mu2 = gv**2 - gu**2 
       guv    = gu*gv
    
       call grid_to_fourier(gv2mu2,fv2mu2,nfx,nfy,ngx,ngy)
       call grid_to_fourier(guv,fuv,nfx,nfy,ngx,ngy)
    
       do jx = 0, nkx
          cjac0(jx,:) = & 
           cmplx(kxky(jx,:)*fv2mu2(2*jx  ,:) + kx2mky2(jx,:)*fuv(2*jx  ,:), &
                 kxky(jx,:)*fv2mu2(2*jx+1,:) + kx2mky2(jx,:)*fuv(2*jx+1,:))
       enddo
          
    end select
    
    if (jac_scale .ne. 1.0 .or. jac_scale .ne. 0.0)  &
                               cjac0(:,:) = jac_scale*cjac0(:,:)
    
    return
    end subroutine jacobian
    
    ! *********************
    ! * SUBROUTINE CQ2FUV *
    ! *********************
    subroutine cq2fuv
    use catmod
    implicit none
    
    integer :: jx
    
    do jx = 0, nkx
       fu(2*jx  ,:) = -aimag(cq(jx,:)) * kyrk2pa(jx,:)
       fu(2*jx+1,:) =  real (cq(jx,:)) * kyrk2pa(jx,:)
       fv(2*jx  ,:) =  aimag(cq(jx,:)) * kxrk2pa(jx,:)
       fv(2*jx+1,:) = -real (cq(jx,:)) * kxrk2pa(jx,:)
    enddo
    
    return
    end subroutine cq2fuv
    
    
    ! ***********************
    ! * SUBROUTINE CQ2FQXQZ *
    ! ***********************
    subroutine cq2fqxqy
    use catmod
    implicit none
    
    integer :: jx, jy
    
    do jx = 0, nkx
       do jy = 0, nfy
          fqx(2*jx  ,jy) =  aimag(cq(jx,jy)) * kx(jx)
          fqx(2*jx+1,jy) = -real (cq(jx,jy)) * kx(jx)
    
          fqy(2*jx  ,jy) =  aimag(cq(jx,jy)) * ky(jy)
          fqy(2*jx+1,jy) = -real (cq(jx,jy)) * ky(jy)
       enddo
    enddo
    
    return
    end subroutine cq2fqxqy
    
    
    ! ***********************
    ! * SUBROUTINE ADD_FORC *
    ! ***********************
    subroutine add_forc
    use catmod
    implicit none
    integer :: i,j,ic,ifk
    real(8) :: k2tmp
    complex(8) :: psif
    real(8) :: eni,enf,ran4
    
    eni=0.d0
    enf=0.d0
    
    select case (nforc)
       case (1) 
          cq = cq + cqfrc
     
       case (2)
          do ifk = 1,nk
             i = in(ifk,1)
             j = in(ifk,2)
    !        psif = ampcoeff*exp(ci*twopi*0.1525125)
             psif = ampcoeff*exp(ci*twopi)
             k2tmp = k2pa(i,j)
             eni = eni+(cq(i,j)*cq(i,j))/k2tmp
             cq(i,j) = cq(i,j)-k2tmp*psif*rnxy
             enf = enf+(cq(i,j)*cq(i,j))/k2tmp
          enddo
    
       case (3)
          ic=tstep+1
          if (mod(ic,itau).eq.0.) then
             call random_number(ran4)
             phi=ran4*twopi*ci
          endif
          psif=ampcoeff*exp(phi)
          ! add the forcing to the spectral vorticity q -> q+ k2tmp *psif
          do ifk=1,nk
             i = in(ifk,1)
             j = in(ifk,2)
             k2tmp = k2pa(i,j)
             eni=eni+(cq(i,j)*cq(i,j))/k2tmp
             cq(i,j) = cq(i,j)-k2tmp*psif*rnxy
             enf=enf+(cq(i,j)*cq(i,j))/k2tmp
          enddo
    
       case (4)
          do ifk=1,nk
             i = in(ifk,1)
             j = in(ifk,2)
             call random_number(ran4)
             psif = ampcoeff*exp(ci*twopi*ran4)
             k2tmp = k2pa(i,j)
             eni=eni+(cq(i,j)*cq(i,j))/k2tmp
             cq(i,j) = cq(i,j)-k2tmp*psif*rnxy
             enf=enf+(cq(i,j)*cq(i,j))/k2tmp
          enddo
    
    end select
    
    return
    end subroutine add_forc
    
    
    ! ************************
    ! * SUBROUTINE ADD_FORC1 *
    ! ************************
    subroutine add_forc1
    use catmod
    implicit none
    
    
    select case (nforc)
       case(1)
    
       case(2)
    
       case(3)
    
       case default
          return
    end select
    
    return
    end subroutine add_forc1
    
    
    
    ! **************************
    ! * SUBROUTINE WRITE_TSERI *
    ! **************************
    subroutine write_tseri
    use catmod
    implicit none
    
    real(8) :: ener,enst,qint
    
    
    qint = sum(gq(:,:))
    ener = 0.5 * rnxy * (sum(gu(:,:) * gu(:,:)) + sum(gv(:,:) * gv(:,:)))
    enst = 0.5 * rnxy *  sum(gq(:,:) * gq(:,:))
    
    write(nutseri,*) tstep,qint,ener,enst
    
    return
    end subroutine write_tseri
    
    
    ! ************************
    ! * SUBROUTINE WRITE_CFL *
    ! ************************
    subroutine write_cfl
    use catmod
    implicit none
    
    real(8) :: cfl
    real(8) :: maxu,maxv
    
    
    !--- check courant number 
    maxu = maxval(abs(gu(:,:))) * dt * ngx / lx
    maxv = maxval(abs(gv(:,:))) * dt * ngy / ly
    
    cfl = max(maxu,maxv)
    
    write(nucfl,*) tstep,cfl
    
    return
    end subroutine write_cfl
    
    
    ! ***********************
    ! * SUBROUTINE INIT_OPS *
    ! ***********************
    subroutine init_ops
    use catmod
    implicit none
    
    integer :: jx,jy    
     
    kx2    = kx*kx
    kx2(0) = 1.0    
    
    ky2    = ky*ky
    
    do jy = 0, nfy
       do jx = 0, nkx
          k2     (jx,jy) = (kx2(jx)+ky2(jy))        ! - Laplacian      psi --> -q
          k2pa   (jx,jy) = (kx2(jx)+ky2(jy))+alpha  ! - mod Laplacian  psi --> -q
          rk2pa  (jx,jy) = 1.0d0/(k2(jx,jy)+alpha)  ! 1/mod Laplacian    q --> -psi
          kxrk2pa(jx,jy) = kx(jx)*rk2pa(jx,jy)      ! q ---> v/i  
          kyrk2pa(jx,jy) = ky(jy)*rk2pa(jx,jy)      ! q ---> u/i
          kx2mky2(jx,jy) = kx2(jx)-ky2(jy)          ! utility Jacobian 3
          kxky   (jx,jy) = kx(jx)*ky(jy)            ! utility Jacobian 3
       enddo
    enddo
    
    return
    end subroutine init_ops
    
    
    ! ***********************
    ! * SUBROUTINE INIT_RST *
    ! ***********************
    subroutine init_rst
    use catmod
    implicit none
    
    inquire(file=cat_rstini,exist=lrst)
    if (lrst) then
      open(nurstini,file=cat_rstini,form='unformatted')
    endif
    
    open(nurstfin,file=cat_rstfin,form='unformatted')
    
    return
    end subroutine init_rst
    
    
    
    ! *************************
    ! * SUBROUTINE CHECK_DIFF *
    ! *************************
    subroutine check_diff(cold,cnew,ytext)
    use catmod
    complex(8) :: cold(0:nkx,0:nfy)
    complex(8) :: cnew(0:nkx,0:nfy)
    character(*) :: ytext
    
    write (nudiag,*) ytext
    write (nudiag,*) "max old  = ",abs(maxval(real(cold)))+abs(maxval(aimag(cold)))
    write (nudiag,*) "max new  = ",abs(maxval(real(cnew)))+abs(maxval(aimag(cnew)))
    write (nudiag,*) "max diff = ",abs(maxval(real (cnew-cold(0:nkx,0:nfy)))) + &
                                   abs(maxval(aimag(cnew-cold(0:nkx,0:nfy))))
    write (nudiag,*) "min old  = ",abs(minval(real(cold)))+abs(minval(aimag(cold)))
    write (nudiag,*) "min new  = ",abs(minval(real(cnew)))+abs(minval(aimag(cnew)))
    write (nudiag,*) "min diff = ",abs(minval(real (cnew-cold(0:nkx,0:nfy)))) + &
                                   abs(minval(aimag(cnew-cold(0:nkx,0:nfy))))
    
    
    return
    end