c %Z% %M% %I% (IDS) %G%

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       this file contains the subroutines corresponding to     c
c       slgcm.main in the 3d model.                             c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      program slgcm
      include "slgcm.slgprm"
      include "slgcm.slgprg"
      include "slgcm.slgcon"
      include "slgcm.mp"

      integer    dt, rdhour
      logical  geopl, dftflg
      common /inout/ dt, initc
CD
      integer icount
      common /serial/ icount
CD

      namelist /inlist/ delt,  mstep, rdhour,  divdmp, ddmpb, ddmpf
     *,                 textt, textn, texpt, texpn, epss, geopl
     *,                 dftflg, krsi, taus, tauc
           
      namelist /outlist/ iohour, krso
      
      common /mean/fmean

      logical first
C
C     Enable nodal program
C
      call mp_enable_nodes()
CD
CD    This variable keeps track of number of interpolations.
CD    This is to ensure that we have a snapshot of the results
CD    after n time steps as well as the last.
      icount = 0
CD
c
c initialize constants
c      
      call consl
      cdf   = 1.0
c
c read namelists
c
      read(51,inlist)
      read(51,outlist)
c
c print namelists
c
      write(6,inlist)
      write(6,outlist)
c
      dt = delt + 0.0001
      cdamp  = divdmp
c
      texf1  = textt
      texf2  = textn
      texfa  = texpt
      texfb  = texpn
c
      dl1    = 0.5 * (1.0 + epss)
      dl2    = 1.0 - dl1
      print *, '    time step will be ',dt,' secs.'

      mmm = 0
      rdt      = dt
      call condt(dt)
c

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     step 1:  specify the initial condition for the shallow water c
c     equation model.                                              c

      if(geopl) then
         print *,'creating artificial initial condition'
         vm = 20.0
         print *,'wind ',vm
         call geopol(vm)
      else
         print *,'reading initial condition from unit ',krsi
         call init(krsi)

      endif
c
c include initial condition in the output file and select output file
c
      ilen=mstep
      call diagout(krso,ilen,'/tmp/2d.out')



      if(dftflg) then
         call dgfiltr(dt,taus,tauc)
         call condt(dt)
         rdt = dt
      endif
      first = .true.

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                  c
c what follows is the actual integration                           c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
100   continue

c     
c     calculate the global kinetic energy (tke), global potential 
c     energy (pot), total energy (eg) and the global mean of 
c     geopotential perturbation (gmass, in meters) which is an
c     indicator of the mass conservation.
c
      call engch(tkeu,tkev,tke,pot,eg,gmass)
      thours = mmm*(dt/abs(dt))

      nmm = mmm*dt/(iohour*3600)
      kmm = nmm*iohour
      lmm = kmm - mmm
      if(lmm.eq.0) then
      mhour = (mmm*dt)/3600

      iii=mhour
      call diagout(krso,iii,'/data1/test1.grads')
      print *,'writing output at ',mmm
      endif
      call step(rdt, first)
      write(6,1111) mmm,tke,pot,eg,gmass
      
      mmm = mmm + 1
      
      if (mmm .ge. mstep) then
         go to 101
      endif

      go to 100

101   continue

CT      call mp_bc_int_from_host(MPQUIT)
CT      call cmmd_sync_host_with_nodes()
CT      call mp_serve_nodes()
CT      call cmmd_sync_host_with_nodes()

      stop         
 1111 format(1x, 'hours=',i4,'  tke=',e11.4, '  pot=',e11.4
     #,          '  eg=',e11.4, '  gmass=',e11.4)

      end
