#include "dims.h"
!
      module init_module
      implicit none
#include "params.h"
!
! Initialize and store model non-input variables and constants:
!
! istep: the current time step index. istep is initialized
! to 0, then incremented before each time step (see advance.f).
!
      integer :: istep     ! time step index
!
! iter (iteration number): the number of timesteps (at the current 
! step length) from model time 0,0,0 to current model time (day,hour,min).
!
      integer :: iter
!
! iyear and iday represent the current model calendar date
! (iyear is 4-digit). Uthr is current decimal hours. These
! are updated per timestep in advnce.
!
      integer :: iyear     ! Current model calendar 4-digit year
      integer :: iday      ! Current model calendar day
      real :: uthr         ! Current ut (decimal hours)
!
! Start date and times are set from input DATE and START time when run 
!   is initial (run is initial when SOURCE and SOURCE_START were provided 
!   at input, but SOURCE_START is ignored wrt start date and time).
! For continuation run, start date and times are read from the startup 
!   history, i.e., they are "carried forward" from the startup history.
!
      integer :: start_year     ! 4-digit starting year 
      integer :: start_day      ! starting day (1-367)
      integer :: start_mtime(3) ! starting model time (day,hr,min)
!
! getgpi flag must go here rather than in gpi_mod to avoid circular 
! module dependency with nchist_mod and init_mod.
!
      integer :: igetgpi   ! 0/1 flag to get GPI data (see gpi_ncfile
                           ! in input_mod.f. Set in init_mod.f)
!
      character(len=16) ::
     |  host,              ! host machine
     |  system,            ! operating system of host (from pre-proc macros)
     |  logname            ! user login name
      character(len=8) ::
     |  rundate,           ! current local date of run
     |  runtime            ! current local time of run
!
! Secs is updated in advnce, and is referenced in opflux, settei,
!   sun, and chapmn.
      real :: 
     |  secs,           ! current ut in seconds
     |  sfeps,          ! solar output change due to orbital eccentricity
!    |  alfalp,efluxlp  ! low-energy protons in aurora
     |  sundec,         ! sun's declination (updated in advnce)
     |  sin_sundec,     ! sin(sundec) (updated in advnceday)
     |  cos_sundec      ! cos(sundec) (updated in advnceday)
!
      contains
!-----------------------------------------------------------------------
      subroutine init
!
! Initialize (this is called by tgcm.F after input):
! Some init also takes place in inp_model (input_mod.f)
!
      use input_module,only: start,step,secflds,secfmag,date,calday,
     |  gpi_ncfile,mxhist_prim,mxhist_sech,output,secout,mkhvols,
     |  nmc,source_start,amienh,amiesh
      use hist_module,only: hist_init,isechist,nsecfmag,nsecfgeo,
     |  nstep,nhist_total,nsech_total,nsource,nseries_prim,nseries_sech,
     |  nfiles_prim,nfiles_sech
      use fields_module,only: set_fprog,set_fsech,set_fsechmag,
     |  fsech,fsechmag
      use cons_module,only: pi,init_cons,set_consdyn,imax
      use amie_module,only: rdamie_nh,rdamie_sh
#include "fgcom.h"
#include "rfft.h"
!
! External:
      integer,external :: mtime_to_nstep
!
! Local:
      real :: theta0
      integer :: i,nfsech
!
! Initialize derived constants (init_cons is in cons_module):
      call init_cons
!
! Initialize dynamo constants 
! (set_consdyn replaces cnstnt.f and is in cons_module)
      call set_consdyn
!
! Initialize buffer indices (init_buff is in this module):
      call init_buff
!
! Get login name:
      logname = ' '
      call getenv('LOGNAME',logname)
      if (len_trim(logname)==0) then
        write(6,"(/,'>>> init: Cannot get LOGNAME environment ',
     |    'variable.',/)")
        stop 'LOGNAME'
      endif
!
! Get host name:
      call gethostsname(host)
!
! Operating system (based on pre-proc macro):
      call setosys(system)
!
! Get run date (current date and time):
      call datetime(rundate,runtime)
!
! Iter is the number of time steps from 0,0,0 to the curent model
! time, using the current step length. Iter is incremented once per 
! timestep in advnce.
!
      iter = mtime_to_nstep(start(:,1),step)
!
! iyear and iday are current calendar year and day.
! If calday > 0, the model is advanced in calendar time, starting
!   at day calday. If calday==0, model is NOT advanced in calendar
!   time and the calendar day is date(2) throughout the run.
! iyear and iday are incremented in advnce if the model is
!   advancing in calendar time.
!
      iyear = date(1)
      iday =  date(2)
      sfeps = 1.
!
! Model is being advanced in calendar time:
! Initialize orbital eccentricity.
      if (calday > 0) then
        iday = calday
        theta0 = 2.*pi*float(iday)/365.
        sfeps = 1.000110+0.034221*cos(theta0)+0.001280*sin(theta0)+
     |          0.000719*cos(2.*theta0)+0.000077*sin(2.*theta0)
      endif
!
! 2/00: these were in modsrc.snoe (tgcm13mt), but were unused.
! Low-energy protons:
!     alfalp = 10.
!     efluxlp = 1.e-20
!
! Set GPI flag:
      igetgpi = 0
      if (len_trim(gpi_ncfile) > 0) igetgpi = 1
      if (igetgpi > 0) 
     |  write(6,"('  gpi_ncfile = ',a)") trim(gpi_ncfile)
!
! ixtimep is 4th dimension index to fg-array for previous time step
! ixtimec is 4th dimension index to fg-array for current time step
! (see fogcm.f)
!
      ixtimep = 1
      ixtimec = 1
!
! Initialize amie, and get amie file if necessary:
!     call init_amie
      if (len_trim(amienh) > 0) then
        write(6,"('Reading  AMIENH = ',a)") trim(amienh)
        call rdamie_nh
      endif
      if (len_trim(amiesh) > 0) then
        write(6,"('Reading  AMIESH = ',a)") trim(amiesh)
        call rdamie_sh
      endif

!
! Initialize non-input history variables for beginning of run:
      call hist_init
!
! setfft calls set99 for fft init. This call returns trigs and ifax
! in global common in rfft.h.
!
      call setfft(trigs,ifax,imax)
!
! Set start date and time if initial run (if continuation run,
! these are set when reading the startup history):
      if (nsource > 0) then ! initial run
        start_year = date(1)
        start_day  = date(2)
        if (calday > 0) start_day = calday
        start_mtime = start(:,1)
      endif
!
! Initialize prognostic fields type:
! (sub set_fprog is in flds_mod.f)
      call set_fprog
!
! Initialize secondary history fields type: 
! Sub set_fsech (flds_mod.f) for fields on geographic grid.
! Sub set_fsechmag (flds_mod.f) for fields on magnetic grid.
! nsecfmag (hist_mod.f), is number of requested magnetic fields,
!   and was set by hist_init, called above.
!
      if (isechist > 0) then
        call set_fsech(secflds,nfsech)
        call set_fsechmag(secfmag)
      endif
!
! Initialize sun's declination:
      sundec=atan(tan(23.5*pi/180.)*sin(2.*pi*float(iday-80)/365.))
      sin_sundec = SIN(sundec) ! C(95)
      cos_sundec = COS(sundec) ! C(96)
!
! Report to stdout:
      write(6,"(/,'Model run initialization:')")
      write(6,"('  nstep      = ',i6,4x,
     |  '(Number of time steps this run)')") nstep
      write(6,"('  iter       = ',i6,4x,
     |  '(Initial iteration number)')") iter
      write(6,"('  iyear      = ',i6,4x,
     |  '(Beginning calendar year)')") iyear
      write(6,"('  iday       = ',i6,4x,
     |  '(Beginning calendar day)')") iday
      write(6,"('  igetgpi = ',i6,4x,
     |  '(If > 0, geophysical indices database will be used.)')")
     |  igetgpi
!
! ncep/nmc are in time-gcm only:
!     write(6,"('  ncep    = ',i6,4x,
!    |  '(If > 0, use NCEP Z and TN 10 mb lower boundaries.')")
!    |  ncep
!     write(6,"('  nmc     = ',i6,4x,
!    |  '(If > 0, use NMC  Z and TN 10 mb lower boundaries.')")
!    |  nmc
!
      if (nsource > 0) then
        write(6,"(/,'This is an initial run:')")
        write(6,"('  start_year = ',i6,5x,
     |    '(Starting year of initial run)')") start_year
        write(6,"('  start_day  = ',i6,5x,
     |    '(Starting day of initial run)')") start_day
        write(6,"('  start_mtime= ',i4,2i3,1x,
     |    '(Starting mtime of initial run)')") start_mtime
      endif
!
! Report re primary histories to stdout:
      write(6,"(/,'Primary Histories:')")
      write(6,"('  nsource       = ',i5,2x,
     |  '(If > 0, a primary source history was provided)')") nsource
      write(6,"('  nseries_prim  = ',i5,2x,
     |  '(Number of primary time series)')") nseries_prim
      write(6,"('  nhist_total   = ',i5,2x,
     |  '(Number of primary histories to be written)')") nhist_total
      write(6,"('  nfiles_prim   = ',i5,2x,
     |  '(Number of primary output files to be written)')") nfiles_prim
      write(6,"('  mxhist_prim   = ',i5,2x,
     |  '(Maximum number of primary histories per file)')") mxhist_prim
!
! Report re secondary histories to stdout:
      if (isechist > 0) then
        write(6,"(/,'Secondary Histories:')")
        write(6,"('  nseries_sech  = ',i5,2x,
     |    '(Number of secondary time series)')") nseries_sech
        write(6,"('  nsech_total   = ',i5,2x,
     |    '(Number of secondary histories to be written)')") nsech_total
        write(6,"('  nfiles_prim   = ',i5,2x,
     |    '(Number of secondary output files to be written)')")
     |    nfiles_sech
        write(6,"('  nsecfmag      = ',i5,2x,
     |    '(Number of secondary history fields on magnetic grid)')") 
     |    nsecfmag
        write(6,"('  mxhist_sech   = ',i5,2x,
     |    '(Maximum number of secondary histories per file)')")
     |    mxhist_sech
!
! Report secondary history fields: 
        write(6,"(/,'Secondary history fields on geographic grid ',
     |    ' (number of fields =',i3,'):')") nsecfgeo
        do i=1,nsecfgeo
          if (fsech(i)%prognostic) then
            write(6,"('  Field ',a,' (',a,')')") fsech(i)%name(1:8),
     |        trim(fsech(i)%long_name)
          else
            write(6,"('  Field ',a,' (diagnostic)')") fsech(i)%name
          endif
        enddo
        write(6,"(/,'Secondary history fields on magnetic ',
     |    'grid (number of fields =',i3,'):')") nsecfmag
        do i=1,nsecfmag
          write(6,"('  Mag field ',a,' (diagnostic)')") fsechmag(i)%name
        enddo
      endif
      end subroutine init
!-----------------------------------------------------------------------
      subroutine init_buff
      use cons_module,only: ncols,kflds,len1,len2,len3,kmax,kmaxp1,
     |  imaxp2
!
! 3/01: Initialize buffer indices, etc. 
! This was in old con.F, and can be eliminated after further rewrite.
!
#include "buff.h"
#include "index.h"
!
! Local:
      integer :: i,k,n,ncbuf,ncphys,nbaphy,lim,maxfld,maxphy,len,
     |  klim
!     DATA KPHYS/nzphys*3/
      integer :: KPHYS(nphys) = (/
     | 3       ,3       ,3         ,3       ,3       ,3        ,
     | 3       ,3       ,3         ,3       ,3       ,3        ,
     | 3       ,3       ,3         ,3       ,3       ,3        ,
     | 3       ,3       ,3         ,3       ,3       ,3        ,
     | 3       ,3       ,3         ,3       ,3       ,3        ,
     | 3       ,3       ,3         ,3       ,3       ,3        ,
     | 3       ,3       ,3/)
!
      LBDSK=0                 ! local (removed from buff.h)
      LBSCM=0                 ! buff.h (needed by rdoldsrc.F)
      LBTAP=1                 ! buff.h (needed by rdoldsrc.F)
      NCOLS=0                 ! ncols is in constant module cons_mod.F
      NCBUF = 0               ! local
      do n=1,nflds                         ! nflds is in params.h
        if (kflds(n)==1) then              ! kflds is in cons_mod
          klim = 1                         ! klim is local
          len = len1                       ! len is local, len1 in cons_mod
        elseif (kflds(n)==2) then
          klim = kmax                      ! kmax in cons_mod
          len = len2                       ! len is local, len2 in cons_mod
        else
          klim = kmaxp1                    ! kmaxp1 in cons_mod
          len = len3                       ! len is local, len3 in cons_mod
        endif
        ndexa(n+1) = ncbuf                 ! ndexa in index.h
        lbscm=lbscm+len
        ncbuf = ncbuf+klim
        if (n.le.ndisk) lbdsk=lbdsk+len    ! ndisk in params.h
        if (n.le.ntape) then
          ncols=ncols+klim
          lbtap=lbtap+imaxp2*klim
        endif
      enddo
      NCPHYS=0                       ! ncphys is local
      do n=1,nphys
        if (kphys(n)==1) then        ! kphys is local
          klim = 1
          len = len1
        elseif (kphys(n)==2) then
          klim = kmax
          len = len2
        else
          klim = kmaxp1
          len = len3
        endif
        ndexb(n+1) = ncphys          ! ndexb in index.h
        ncphys = ncphys+klim
      enddo
C     ****                   COMPUTE COLUMN POINTERS FOR SCM BUFFERS
      nbaddr(1) = 1          ! nbaddr in buff.h, used by init and advnce
      do i=2,8
        nbaddr(i) = nbaddr(i-1)+ncbuf
      enddo
      nbaphy = nbaddr(8)+ncbuf ! nbaphy is local 
      do n=1,nphys
        ndexb(n+1) = ndexb(n+1)+nbaphy
      enddo
C     ****                   CHECK IF F DIMENSION IS APPROPRIATE
      MAXFLD=ZFLDX           ! local
      MAXPHY=ZPHYX           ! local
      write(6,"(/,'F-array dimensioned for fldx=',i5,' phyx=',i5,
     | /,' Model needs ncbuf=',i5,' ncphy=',i5,/)")
     |  maxfld,maxphy,ncbuf,ncphys
      if (ncbuf > maxfld .or. ncphys > maxphy) then
        write(6,"(/,'Exit called by init_buff: need fldx.ge.',i5,
     |    ' phyx.ge.',i5)") ncbuf,ncphys
        stop 'init_buff'
      endif
      end subroutine init_buff
      end module init_module