!=============================================================================
!
! $Id: uci_control.F90,v 1.2 2007/05/11 01:29:39 kouatch Exp $
!
! CODE DEVELOPER
!   Deanne Proctor, et. al., LLNL (Original code from Michael Prather, UCI)
!   ddproctor@llnl.gov
!
! FILE
!   uci_control.F
!
! ROUTINES
!   Uci_Control
!   Set_Vert_Uci
!   Wind
!   Exitc
!   Jshift1
!   Ad_To_Mass
!   Jshiftn
!   Stt_To_Const
!   Get_Diag_1d_Dump_Flag
!
!=============================================================================


!-----------------------------------------------------------------------------
!
! ROUTINE
!   Uci_Control
!
! DESCRIPTION
!   This routine is the wrapper for UCI 3D transport model.
!
!   This transport driver is based on the GISS GCM.
!   This version is meant to be used for both trop and strat:
!     GISS 4x5 / 8x10 x 9 layer GCM history tapes with 5-day averages,
!     GISS 7.8x10 x 23    layer GCM history tapes with monthly convection.
!
! ARGUMENTS
!   const_ncar_infile  : const input file name
!   met_uci_infile     : met   input file name (UCI format)
!c  pr_diag_1d         : should UCI 1D diagnostic output be written?
!c  pr_diag_2d         : should UCI 2D diagnostic output be written?
!c  time_for_diag_init : time for UCI diagnostics initialization?
!   gmi_lm             : UCI vertical grid dimension
!   num_time_steps     : number of time steps
!   nymd               : year/month/day (YYYYMMDD)
!   nhms               : hour/min/sec   (HHMMSS)
!   gmi_sec            : total Gmimod seconds (s)
!   tdt                : model time step (s)
!   uci_lat_out        : UCI latitudes (deg)
!   press              : midpoint  pressures by level (mb)
!   prslay             : interface pressures by level (mb)
!   mass               : total mass of the atmosphere within each grid box
!                        (kg)
!   const              : species concentration, known at zone centers
!                        (mixing ratio)
!
!-----------------------------------------------------------------------------


      subroutine Uci_Control  &,53
!c   &  (const_uci_infile, met_uci_infile, pr_diag_1d, pr_diag_2d,
!c   &   time_for_diag_init, gmi_lm, num_time_steps, nymd, nhms,
     &  (const_uci_infile, met_uci_infile, gmi_lm, num_time_steps,  &
     &   nymd, nhms, gmi_sec, tdt, uci_lat_out, press, prslay,  &
     &   mass, const)

!xx!? implicit none

#     include "netcdf.inc"
#     include "gmi_dims.h"
#     include "uci.h"


!     ----------------------
!     Argument declarations.
!     ----------------------

      character*(*) :: const_uci_infile
      character*(*) :: met_uci_infile
!c    logical :: pr_diag_1d, pr_diag_2d
!c    logical :: time_for_diag_init
      integer :: gmi_lm
      integer :: num_time_steps
      integer :: nymd, nhms
      real*8  :: gmi_sec
      real*8  :: tdt
      real*8  :: uci_lat_out(ju1_gl:j2_gl)
      real*8  :: press (k1:k2)
      real*8  :: prslay(k1-1:k2)
      real*8  :: mass  (i1:i2, ju1:j2, k1:k2)
      real*8  :: const (i1:i2, ju1:j2, k1:k2, num_species)


!     ----------------------
!     Variable declarations.
!     ----------------------

      character*78  :: titver

      character*128 :: uci_read_file
      character*128 :: uci_write_file

      character*8   :: tlabl(10)

      logical       :: dump_diag
      logical, save :: first_time = .true.

      integer :: cdfid
      integer :: dnum_time_steps
      integer :: dnymd, dnhms
!xx   integer :: idt
      integer :: ierr_getlun
      integer :: ijx
      integer :: il, ij, ik, ic
      integer :: gmi_rec

      real*8  :: gmi_day

      real*8  :: mw(num_species)

      real*8  :: ad_jsh   (i1:i2, ju1:j2, k1:k2)
      real*8  :: const_jsh(i1:i2, ju1:j2, k1:k2, num_species)

!     ====
      save
!     ====

      data  tlabl / 'INITIAL ', 'WINDS+P ', 'DYNAMICS', 'ALL-CONV',  &
     &              '________', 'SOURCES ', 'CHEMSTRY', 'DIFF-MIX',  &
     &              'DIAGNOST', 'SAVETAPE' /

      data  titver  &
     &  / ' CTM (GISS/Harvard/UCIrvine) p-1.0 1/95(Prather)' /


!     ----------------
!     Begin execution.
!     ----------------

      if (pr_diag) then

        if (first_time) then

          Write (6,*) 'Uci_Control called for the first time.'

        else

          Write (6,*) 'Uci_Control called, ', ntau * 3600.0d0

        end if

      end if


      if (.not. first_time) then
!       =======
        goto 20
!       =======
      end if


!rt
!rt
!rt
!rt          first time through UCI
!rt
!rt
!rt
!xx   ymd=start_ymd
!xx   hms=start_hms
!xx
      KYY = 93
      KMM = 7
      KDD = 1
!     CALL IDATE(KMM,KDD,KYY)
!     CALL SYSTEM('date')
!---------------------input files--(___.ctm)----------------------------
!#1         CTM INITIALIZE/RESTART FILE
!xx   OPEN(1,FILE='init.uci',status='OLD',form='UNFORMATTED')
!#2         ALTITUDE/LAND-FRACT FOR GCM GRID
!xx   OPEN(2,FILE='landalt.uci',STATUS='OLD',form='FORMATTED')
!#3         GCM DATA TAPE (winds, etc.)
!rt
!rt   OPEN(3,FILE='wind.ctm',STATUS='OLD',FORM='UNFORMATTED')
      call opencdf(met_uci_infile,cdfid)
!rt
!rt
!xx
      call FI_Getlun (lun8, ierr_getlun)
!xx   OPEN(8,FILE='input.uci',STATUS='OLD',form='FORMATTED')
      OPEN(unit=lun8,FILE='input.uci',STATUS='OLD',form='FORMATTED')
!xx
!#9         TRACER SETUP DATA CONC/SOURCE/CHEM
      call FI_Getlun (lun9, ierr_getlun)
!xx   OPEN(9,FILE=const_uci_infile,STATUS='OLD',form='FORMATTED')
      OPEN(unit=lun9,FILE=const_uci_infile,  &
     &     STATUS='OLD',form='FORMATTED')
!--------------------output files--(ctm.___)----------------------------
!#10        (FORMATTED DIAGNOSTIC OUTPUT, 80 col)
      call FI_Getlun (lun10, ierr_getlun)
!xx   OPEN(10,FILE='uci.pch',form='FORMATTED')
      OPEN(unit=lun10,FILE='uci.pch',form='FORMATTED')
!#11        (CTM DAILY RESTART FILE, UNFORMATTED)
!xx   OPEN(11,FILE='uci.sav',form='UNFORMATTED')
!#12        (UNFORMATTED DIAGNOSTICS)
      call FI_Getlun (lun12, ierr_getlun)
!xx   OPEN(12,FILE='uci.unf',form='UNFORMATTED')
      OPEN(unit=lun12,FILE='uci.unf',form='UNFORMATTED')
!xx
      call FI_Getlun (lun13, ierr_getlun)
!xx   OPEN(13,FILE='output.uci',form='FORMATTED')
      OPEN(unit=lun13,FILE='output.uci',form='FORMATTED')
!xx
!-----------------------------------------------------------------------
!xx
      WRITE (lun13,909)
      WRITE (lun13,900) TITVER,TITCOM,KDD,KMM,KYY
      WRITE (lun13,909)
!xx
!rt
!rt
!rt
!rt          name UCI NetCDF restart file
!rt
      UCI_read_file  = "UCI_read"
      UCI_write_file = "UCI_write"
!rt
!rt
!rt
!xx   CALL INPUT
      CALL INPUT (tdt,UCI_read_file)

!xx
      ijx = 1

      do ij = ju1_gl, j2_gl
        uci_lat_out(ij) = uci_lat(ijx)
        ijx = ijx + 1
      end do
!xx

!rt
!xx
!
      CALL INPTR
!
!
      CALL GRDSET
      WRITE(lun10,981) (XLABEL(K),K=1,15),NTRACE,IM,JM,LM,I0,J0,IYEAR
      WRITE(lun12)     (XLABEL(K),K=1,15),NTRACE,IM,JM,LM,I0,J0,IYEAR
!xx   CLOSE (1)
!xx   CLOSE (2)
      CLOSE (lun9)
      KDACC = 0
      KWACC = 0
      DO K=1,KACCZ
        SAVACC(K) = 0.0
      ENDDO
!-----------------------------------------------------------------------
!-------Begin CTM calculation at time  = NTAU hrs relative to IYEAR-----
!-------NTAU = ABSOLUTE TIME (I*4), ITAU = RELATIVE TIME (0 - <8760 HRS)
!-------TAU (R*4) = NTAU, calculation not resolved finer than 1 hour
      NTAU   = (TAU+0.1)
      NTAUE  = (TAUE+0.1)
      IF (NTAU .GE. NTAUE) GOTO 99
      ITAU   = MOD(NTAU,8760)
!-----------------------------------------------------------------------
!xx   CALL WIND (0,ITAU,ITAUW)
      CALL WIND (0,ITAU,ITAUW,cdfid)
!-----------------------------------------------------------------------
      IF (NDYN*(NREAD/NDYN) .NE. NREAD ) GOTO 97
      NTAU   = (NTAU/NREAD)*NREAD
      CALL DAILY
      LPRT   = ND70.GT.0
      WRITE (lun13,901) TLABL(1),NTAU,JDATE,JMONTH,JYEAR
      NTAU00 = NTAU
!rt
      IF (LCONT)  THEN
        NTAU00  = TAU0
        TAU0    = TAU
      ENDIF
!
      if (LCONT)      Write (lun13,905) NTAU,NTAU00,TAU0,TAU
      if (.not.LCONT) Write (lun13,906) NTAU,NTAU00,tAU0,TAU
!rt
!
      NTSAVE = -9999
!-----------------------------------------------------------------------
!
      CALL AIRSET
!
      CALL TRACE0
!
!-----------------------------------------------------------------------
!rt
!rt      IF (LCONT) THEN
!rt        CALL DIAG0(9)
!rt      ELSE
           CALL DIAG0(1)

!          should really pass mw in from Gmimod?  !!
           do k=1,num_species
             mw(k)=28.97/tcvv(k)
           end do
!xx        call Get_Diag_1d_Prnt_Flag (pr_diag_1d)
!xx        call Get_Diag_2d_Prnt_Flag (pr_diag_2d)

!c         if(pr_diag_1D) then
!             need j shifts in const and ad being fed to diag0_gmi
!c            call jshiftn(stt,const_jsh)
!c            call jshift1(ad,ad_jsh)
!c            call diag0_gmi (-1,const_jsh,ad_jsh)
!c            call jshiftn(stt,const_jsh)
!c            call jshift1(ad,ad_jsh)
!c            call diag0_gmi (1,const_jsh,ad_jsh)
!c         endif

!c         if(pr_diag_2d) then
!             note there's no shift for Gmi_Diag_Save so should
!             rename const_jsh,ad_jsh when get rid of 1d diagnostics
!c            call stt_to_const(stt,const_jsh)
!c            call ad_to_mass(ad,ad_jsh)
!c            call Gmi_Diag_Save
!c   &          (1,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c         endif

!rt      ENDIF
!rt
      MONTH0 = 0
!rt
!rt
!rt
!                  set month_co2
!                  to first of year relative to IYEAR input value
!                  month_co2 = 0, 12, 24, ... ,180
!
      if (run_type .eq. " CO2") then
      month_co2 = (JYEAR - IYEAR)*12
      endif
!rt
!rt
!rt
      DTWIND = 3600.D0 *REAL(NDYN)
      DTCONV = 3600.D0 *REAL(NCONV)

!xx
!     =================
      call Set_Vert_Uci  &
!     =================
     &  (k1, k2, gmi_lm, ltm, ptop, sige, plev, press, prslay)

      first_time = .false.
!xx


!-----------------------------------------------------------------------
!------------------------MAIN LOOP--------------------------------------
!------------increment by NDYN hours as basic step----------------------
!-----------------------------------------------------------------------

!     ========
 20   continue
!     ========

!rt
!rt
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,914) NDYN,ITAU,ITAUW,NTAU
      endif
!rt
!rt
!rt
!
!---------do not read new winds if NDYN < NREAD and keep using same winds
      IF (ITAU+NDYN.GT.ITAUW) THEN
        CALL WIND (1,ITAU,ITAUW,cdfid)
        CALL DYN0(1)

!       in principle only have to fill mass array when you want to
!       have pr_ascii1_c true, but...
!       (also assumes only change mass when get new windfield+dyn0
        call jshift1(ad,mass)
!       could be using mass below instead of ad_jsh?

        CALL DIAG0(9)
!c      if(pr_diag_1d) then
!c        call jshiftn(stt,const_jsh)
!c        call jshift1(ad,ad_jsh)
!c        call diag0_gmi (9,const_jsh,ad_jsh)
!c      endif

        IF (LPRT) WRITE(lun13,901) TLABL(1),NTAU,JDATE,JMONTH,JYEAR
      ENDIF
!----------CHECK FOR DAILY / MONTHLY HOUSEKEEPING
      IF (MOD(NTAU,24).EQ.0) THEN
        CALL DAILY
      ENDIF
      IF (MONTH .NE. MONTH0) THEN
        MONTH0 = MONTH
!rt
!rt
!rt
      if ((run_type .eq. " N2O") .or.  &
     &    (run_type .eq. " F11")) then
!rt
        CALL STRATL
!rt
      endif
!rt
!rt
!rt
      if (run_type .eq. " CO2") then
        month_co2 = month_co2 + 1
!rt
        if (edit_type .eq. "debug") then
          Write (lun13,916) month0,month_co2,tparm_co2(1,10,month_co2,1)
        endif
      endif
!rt
!rt
!rt
      ENDIF
!rt
!rt
!rt
      if (run_type .eq. " CO2") then
!rt
        CALL STRATL_co2
!rt
      endif
!rt
!rt
!rt
      if (run_type .eq. " SF6") then
!rt
        CALL STRATL_sf6
!rt
      endif
!rt
!rt
!rt
!-----------------------------------------------------------------------
!-----Begin sequence of operations on the tracer distribution:
!-----  currently the order is SOURCE/CHEM/ADVECT/PFILT/DIFF/CONV/DIAG
!-----  N.B. the order may affect diagnostics, esp. for "rapid change"
!-----       species with large source/sink terms.
!-----This loop calculates changes in tracer from NTAU to NTAU+NDYN
!-----------------------------------------------------------------------
!xx       call Get_Diag_Init_Flag(time_for_diag_init)
!c        if(time_for_diag_init) then
!c           call stt_to_const(stt,const_jsh)
!c           call ad_to_mass(ad,ad_jsh)
!c           call Gmi_Diag_Save
!c   &       (1,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c        endif
               CALL DIAG0(9)
!c             if(pr_diag_1d) then
!c                call jshiftn(stt,const_jsh)
!c                call jshift1(ad,ad_jsh)
!c                call diag0_gmi (9,const_jsh,ad_jsh)
!c             endif
!----------SOURCE TERMS:------------------------------------------------
      IF (MOD(NTAU,NSRCE).EQ.0) THEN
!rt
!rt
!rt
      if ((run_type .eq. " N2O") .or.  &
     &    (run_type .eq. " F11")) then
!rt
        CALL SOURCE
!rt
      endif
!rt
!rt
!rt
        CALL DIAG0(6)
!c      if(pr_diag_1d) then
!c         call jshiftn(stt,const_jsh)
!c         call jshift1(ad,ad_jsh)
!c         call diag0_gmi (3,const_jsh,ad_jsh)
!c      endif
!c      if(pr_diag_2d) then
!c         call stt_to_const(stt,const_jsh)
!c         call ad_to_mass(ad,ad_jsh)
!c         call Gmi_Diag_Save
!c   &      (6,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c      endif


        IF (LPRT) WRITE (lun13,901) TLABL(6),NTAU,JDATE,JMONTH,JYEAR
      ENDIF
!----------CHEMISTRY:---------------------------------------------------
      IF (MOD(NTAU,NCHEM).EQ.0) THEN
!
        CALL CHEM
!
        CALL DIAG0(7)
!c      if(pr_diag_1d) then
!c         call jshiftn(stt,const_jsh)
!c         call jshift1(ad,ad_jsh)
!c         call diag0_gmi (4,const_jsh,ad_jsh)
!c      endif
!c      if(pr_diag_2d) then
!c         call stt_to_const(stt,const_jsh)
!c         call ad_to_mass(ad,ad_jsh)
!c         call Gmi_Diag_Save
!c   &      (3,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c      endif

        IF (LPRT) WRITE (lun13,901) TLABL(7),NTAU,JDATE,JMONTH,JYEAR
      ENDIF
!----------ADVECT TRACERS BY THE WINDS:  U-V-W pattern based on 8x10 winds
      IF (MOD(NTAU,NDYN).EQ.0) THEN
        DTALFA = DTWIND/4.
        DTBETA = DTWIND/3.
        DTGAMA = DTWIND/2.
!
        CALL DYN2U(DTALFA)
!
        CALL DYN2V(DTBETA)
!
        CALL DYN2W(DTGAMA)
!
        CALL DYN2U(DTALFA)
!
        CALL DYN2V(DTBETA)
!
        CALL DYN2U(DTALFA)
!
        CALL DYN2W(DTGAMA)
!
        CALL DYN2V(DTBETA)
!
        CALL DYN2U(DTALFA)
!
        CALL DIAG0(3)
!c      if(pr_diag_1d) then
!c         call jshiftn(stt,const_jsh)
!c         call jshift1(ad,ad_jsh)
!c         call diag0_gmi (5,const_jsh,ad_jsh)
!c      endif
!c      if(pr_diag_2d) then
!c         call stt_to_const(stt,const_jsh)
!c         call ad_to_mass(ad,ad_jsh)
!c         call Gmi_Diag_Save
!c   &      (5,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh,)
!c      endif
!
!
!       WRITE(lun13,'(11HDYN2. LD42=,I3)') LD42
!       DO J = 14,10,-1
!         WRITE(lun13,'(1P10E8.1)') (AIJ(I,J,1+ND42),I=1,10)
!       ENDDO
        IF (LPRT) WRITE (lun13,901) TLABL(3),NTAU,JDATE,JMONTH,JYEAR
      ENDIF
!----------DIFFUSION OF TRACER:-----------------------------------------
      IF (MOD(NTAU,NDYN).EQ.0) THEN
        CALL DIFFU
        CALL DIFFV
        CALL DIFFW
        CALL DIAG0(8)
!c      if(pr_diag_1d) then
!c         call jshiftn(stt,const_jsh)
!c         call jshift1(ad,ad_jsh)
!c         call diag0_gmi (6,const_jsh,ad_jsh)
!c      endif
!c      if(pr_diag_2d) then
!c         call stt_to_const(stt,const_jsh)
!c         call ad_to_mass(ad,ad_jsh)
!c         call Gmi_Diag_Save
!c   &      (7,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c      endif
        IF (LPRT) WRITE (lun13,901) TLABL(8),NTAU,JDATE,JMONTH,JYEAR
      ENDIF
!----------CONVECTION:  includes boundary layer, moist cumulus, large-scale
      IF (MOD(NTAU,NCONV).EQ.0) THEN
        CALL CNVWET (DTCONV)
        CALL DIAG0(4)
!c      if(pr_diag_1d) then
!c         call jshiftn(stt,const_jsh)
!c         call jshift1(ad,ad_jsh)
!c         call diag0_gmi (7,const_jsh,ad_jsh)
!c      endif
!c      if(pr_diag_2d) then
!c         call stt_to_const(stt,const_jsh)
!c         call ad_to_mass(ad,ad_jsh)
!c         call Gmi_Diag_Save
!c   &      (8,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c      endif
        CALL CNVDRY (DTCONV)
        CALL DIAG0(5)
!c      if(pr_diag_1d) then
!c         call jshiftn(stt,const_jsh)
!c         call jshift1(ad,ad_jsh)
!c         call diag0_gmi (2,const_jsh,ad_jsh)
!c      endif
!c      if(pr_diag_2d) then
!c         call stt_to_const(stt,const_jsh)
!c         call ad_to_mass(ad,ad_jsh)
!c         call Gmi_Diag_Save
!c   &      (8,nymd,nhms,gmi_sec,mw,ad_jsh,const_jsh)
!c      endif
        IF (LPRT) WRITE (lun13,901) TLABL(4),NTAU,JDATE,JMONTH,JYEAR
!       WRITE(lun13,'(11HCONV. LD42=,I3)') LD42
!       DO J = 14,10,-1
!         WRITE(lun13,'(1P10E8.1)') (AIJ(I,J,1+ND42),I=1,10)
!       ENDDO
      ENDIF
!-----------------------------------------------------------------------
!----------COMPLETED A TIME STEP:  update NTAU by NDYN (basic time step)
!-----------------------------------------------------------------------
      NTAU   = NTAU+NDYN
!rt
!rt
!rt
      ITAU   = MOD(NTAU-1,8760)+1
!rt
!rt
!rt
      TAU    = NTAU
      IDAY   = 1+NTAU/24
      TOFDAY = TAU-(IDAY-1)*24
!----------DIAGNOSTICS:  ACCUMULATE SOME DIAGNOSTICS AT END OF TIME STEP
      IF (MOD(NTAU,NDIAG).EQ.0) THEN
        CALL DIAG1
        IF (LPRT) WRITE (lun13,902) NTAU
      ENDIF
      IF (MOD(NTAU,NINST).EQ.0) CALL DIAG6
!----------DIAGNOSTICS: print / zero-out accumulating arrays:-----------
      IF (MOD(NTAU,24).EQ.0) THEN
        CALL DAILY
!       njday is table (input) with value 3 on days when dump wanted
        JDO = NJDAY(JDAY)
        LZERO = MOD(JDO,2)
        LPRNT = MOD(JDO/2,2)
!----------WRITE DIAGNOSTICS to UNIT=6 (printout) and UNIT=10 (pch)-----
        IF (LPRNT.EQ.1) THEN
          CALL DIAG0(9)
!c        if(pr_diag_1d) then
!c           call jshiftn(stt,const_jsh)
!c           call jshift1(ad,ad_jsh)
!c           call diag0_gmi (9,const_jsh,ad_jsh)
!c        endif
          CALL DIAG0(99)
          CALL DIAG2
          CALL DIAG4
          CALL DIAG3
          CALL DIAG5
!---------may need:  ENDFILE 10 / BACKSPACE 10 to clear buffer if crash
        ENDIF
!---------ZERO OUT AND RESET ACCUMLATING ARRAYS-------------------------
        IF (LZERO.EQ.1) THEN
          TAU0  = TAU
          NTAU0 = NTAU
          IDAY0 = IDAY
          TOFDY0= TOFDAY
          JDATE0= JDATE
          JMNTH0= JMONTH
          JYEAR0= JYEAR
          KDACC = 0
          KWACC = 0
          DO 50 K=1,KACCZ
   50     SAVACC(K) = 0.0
          CALL DIAG0(1)
        ENDIF
      ENDIF

!     GMIMOD DIAGNOSTICS DUMP?

!c    if(pr_diag_1d) then
!        since Gmimod ticks clock then checks time, in order to sync
!        diagnostics need to get dump flag for incremented ymd
!c       dnum_time_steps=num_time_steps
!c       dgmi_sec=gmi_sec
!c       dnymd=nymd
!c       dnhms=nhms
!xx      idt=Nint(tdt)
!xx      call Gmi_Tick(dnum_time_steps,dgmi_sec,dnymd,dnhms,idt,tdt)
!c       call Gmi_Tick
!c   &     (dnum_time_steps,dgmi_sec,dnymd,dnhms,tdt)
!c       call Get_Diag_1d_Dump_Flag (dnymd,dump_diag)
!c       if(dump_diag) then
!c           call jshiftn(stt,const_jsh)
!c           call jshift1(ad,ad_jsh)
!c           call diag0_gmi (99,const_jsh,ad_jsh)
!c           call diag0_gmi (-1,const_jsh,ad_jsh)
!c           call diag0_gmi ( 1,const_jsh,ad_jsh)
!c       endif
!c    endif

!      END GMIMOD DIAGNOSTICS DUMP


!----------SAVE COMPLETE CTM CALCULATION: for restart or re-init.-------
!rt   IF (MOD(NTAU,NWRITE).EQ.0) THEN
!xx     REWIND 11
!xx     WRITE(11) TAU,STT,SUT,SVT,SWT,AD,SAVLAB,SAVDAT,SAVACC,SAVTRA
!xx     IF (LSOM) WRITE(11) TAU,SUU,SVV,SWW,SUV,SUW,SVW
!rt
!rt
!rt
!rt                Write NetCDF restart file every NWRITE hours
!rt
      if (mod(NTAU,NWRITE).eq.0) then
        gmi_rec=1
        gmi_day=NTAU/24
        N=1
!rt
        call Write_Mpra2 (UCI_write_file, gmi_rec, gmi_day,  &
     &  STT, ad, SUT, SVT, SWT, SUU, SVV, SWW, SUV, SUW, SVW,  &
     &  IM, JM, LM, N,  &
     &  NTAU00, TAU)
!rt
      if (LCONT)      Write (lun13,907) NTAU,NTAU00,TAU0,TAU
      if (.not.LCONT) Write (lun13,908) NTAU,NTAU00,tAU0,TAU
!rt
!rt
!rt
        NTSAVE = NTAU
        IF (LPRT) WRITE (lun13,901) TLABL(10),NTSAVE,JDATE,JMONTH,JYEAR
!xx     REWIND 11
      ENDIF
!-----------------------------------------------------------------------
!----------TEST FOR TERMINATION OF RUN
!-----------------------------------------------------------------------
!xx
!xx   IF (NTAU.LT.NTAUE) GOTO 20
!     ========
      goto 888
!     ========
!xx
!-----------------------------------------------------------------------
!----------END OF MAIN LOOP, CLOSE DOWN AND SAVE RESULTS----------------
!-----------------------------------------------------------------------
!xx   IF (NTAU.NE.NTSAVE) THEN
!xx     Write(lun13,*) ' begin to write ctm.sav file'
!xx     REWIND 11
!xx     WRITE(11) TAU,STT,SUT,SVT,SWT,AD,SAVLAB,SAVDAT,SAVACC,SAVTRA
!xx     IF (LSOM) WRITE(11) TAU,SUU,SVV,SWW,SUV,SUW,SVW
!xx     REWIND 11
!xx     WRITE (lun13,901) TLABL(10),NTAU,JDATE,JMONTH,JYEAR
!xx   ENDIF
!xx   GOTO 99
!-----------------------------------------------------------------------
!---------PROBLEM WITH MONTHLY FILE-------------------------------------
!xx96 WRITE (lun13,961) NTAU,IMON
!xx   WRITE (lun13,961) NTAU
!xx   GOTO 99
!---------INVALID TIME STEPS--------------------------------------------
   97 WRITE(lun13,971) NDYN,NREAD
      GOTO 99
!-------------------------FINAL-----------------------------------------
   99 CONTINUE
      WRITE (lun13,909)
      WRITE (lun13,900) TITVER,TITCOM,KDD,KMM,KYY
      WRITE (lun13,910) NTAU00,NTAU,JDATE,JMONTH,JYEAR
!     CALL SYSTEM('date')
      WRITE (lun13,909)
      STOP
!xx
 888  continue
!rt
!rt
!rt
      do 890 ic = 1, num_species
        do 890 ik = k1, k2
          do 890 ij = ju1, j2
            do 890 il = i1, i2
              const(il,ij,ik,ic) = stt(il,ij+1,ik,ic)
              const(il,ij,ik,ic) = TCVV(ic) * const(il,ij,ik,ic)
              const(il,ij,ik,ic) = const(il,ij,ik,ic)/ad(il,ij+1,ik)
!rt
 890  continue
!rt
!rt
!rt
      return
!xx
!-----------------------------------------------------------------------
  900 FORMAT(1X,A80/'      COMMON BLOCKS:  ',A40,'    dd/mm/yy=',3I3)
  901 FORMAT(1X,A8,'  TAU=',I10,I6,A4,I5)
  902 FORMAT(' DIAGNOSTICS ACCUMULATED  TAU=',I10,10I4)
!c904 FORMAT('  READ 5-DAY/MONTHLY STATS:  ',3I10)
  909 FORMAT(1X,32('****'))
  910 FORMAT (' CTM finish: NTAU0/NTAU=',2I10,I6,A5,I5)
!c911 FORMAT(' END OF GCM WINDS AT TAUW=',F10.1)
!c912 FORMAT(' REWOUND GCM WIND FILE',2I10)
!c921 FORMAT(' INCORRECT FILE???  NTAU,NTAUW,NTAUC,NREAD=',4I15)
!c931 FORMAT(' =====ERROR ON WIND/CTM FILES')
!c951 FORMAT(' =====UNEXPECTED EOF ON WINDS',F10.1)
!c961 FORMAT(' =====ERROR ON MONTHLY FILES',2I10)
  961 FORMAT(' =====ERROR ON MONTHLY FILES',I10)
  971 FORMAT(' =====ERROR ON NDYN/NREAD:',2I10)
  981 FORMAT(15A4,6I2,I8)
!rt
!rt
!rt
 905  format (1x,"restart",2x,"NTAU=",i8,1x,"NTAU00=",i8,  &
     &   1x,"TAU0=",e12.5,1x,"TAU=",e12.5)
 906  format (1x,"  start",2x,"NTAU=",i8,1x,"NTAU00=",i8,  &
     &   1x,"TAU0=",e12.5,1x,"TAU=",e12.5)
 907  format (1x,"restart (write)",2x,"NTAU=",i8,1x,"NTAU00=",i8,  &
     &   1x,"TAU0=",e12.5,1x,"TAU=",e12.5)
 908  format (1x,"  start (write)",2x,"NTAU=",i8,1x,"NTAU00=",i8,  &
     &   1x,"TAU0=",e12.5,1x,"TAU=",e12.5)
 914  format (1x,"MPRA - CONTROL",1x,"NDYN=",i2,1x,"ITAU=",i8,  &
     &   1x,"ITAUW=",i8,1x,"NTAU=",i8)
 916  format (1x,"CALL STRATL",1x,"month0=",i3,1x,"month_co2=",i3,  &
     &   1x,"tparm_co2=",e12.5)
!rt
!rt
!rt
      END


!-----------------------------------------------------------------------------
!
! ROUTINE
!   Set_Uci_Vert
!
! DESCRIPTION
!   This routine sets Gmimod press & prslay for 2D diagnostics (this is an
!   extra layer of stuff since 2D diagnostics were originally written for
!   UCI, but....).
!
! ARGUMENTS
!   k1, k2   : Gmimod vertical grid range
!   gmi_lm   : UCI vertical grid dimension
!   ltm      : top level still defined by sigma coordinates at lower boundary
!   ptop     : pressure at trop/strat interface (mb)
!   sige     : UCI sigma,
!                sige = (p-ptop)/(psf-ptop) in strat,
!                sige = p/ptop in trop?
!   plev_uci : midpoint  pressures by level (1) &
!              interface pressures by level (2) (mb)
!   press    : midpoint  pressures by level (mb)
!   prslay   : interface pressures by level (mb)
!
!-----------------------------------------------------------------------------


      subroutine Set_Vert_Uci  & 1
     &  (k1, k2, gmi_lm, ltm, ptop, sige, plev_uci, press, prslay)

      implicit none

!c#     include "gmi_diag_constants_uci.h"


!     ----------------------
!     Argument declarations.
!     ----------------------

      integer  k1, k2
      integer  gmi_lm
      integer  ltm

!c?   For tmp fix.
      real*8   ptop
      real*8   sige(*)

      real*8   plev_uci(22,2)

      real*8   press (k1:k2)
      real*8   prslay(k1-1:k2)


!     ----------------------
!     Variable declarations.
!     ----------------------

      integer ik


!     ----------------
!     Begin execution.
!     ----------------

      press(k1:k2) = plev_uci(1:k2,1)


      do ik = 1, ltm+1
        prslay(ik-1) = sige(ik) * (PSF - ptop) + ptop
      end do

      do ik = ltm+2, gmi_lm+1
        prslay(ik-1) = sige(ik) * ptop
      end do


      return

      end


!xx   SUBROUTINE WIND(ISW,ITAU,ITAUW)

      SUBROUTINE WIND(ISW,ITAU,ITAUW,cdfid) 2,9
!
!----ISW=0  = locate pressure fields for initialization of dry-air mass
!----ISW>0  = read in next wind, pressure, temperature, q fields.
!----ITAU   = starting time (hours, 0-8759) of initial calculation,
!----         used to position file   (ONLY used for init: ISW=0)
!----ITAUW  = hours (1-8760) of last read wind record, i.e.,
!----         pressure at ITAUW, winds & q averaged over ITAUW-NREAD to ITAUW
!-----------------------------------------------------------------------
#        include "uci.h"
!xx
      integer cdfid
!xx
!rt
!rt
!rt
      integer ihyper
      integer number
!rt
!rt
!rt
!xx$$
      save
!xx$$
!-----------------------------------------------------------------------
      IF (ISW .GT. 0) GOTO 41
!------------locate first wind/press/etc field--------------------------
!         NTBEG/NTEND = START/END HOURS OF WIND TAPE
!         NAVG        = # HOURS FOR AVERAGE STATS ***NOT USED NOW***
!         NREAD       = # HOURS FOR WINDS (4 OR 8 HOURS USUALLY)
!         II0/IIM     = GRID OFFSET/SIZE IN THE W -> E DIRECTION
!         JJ0/JJM     =                         S -> N
!         LL0/LLM     =                        VERTICAL
!*********THE OPTION (II0,JJ0,LL0) IS NO LONGER ALLOWED,****************
!*********MUST HAVE GLOBAL GRID OF WINDS:  II0=JJ0=LL0=0****************
!-----------------------------------------------------------------------
!rt   READ(3) NTBEG,NTEND,NAVG,NREAD,I00,IIM,J00,JJM,L00,LLM
!rt
!rt
!rt
      call met_readh (cdfid,number,  &
     &   NTBEG,NTEND,NAVG,NREAD,  &
     &   I00,IIM,J00,JJM,L00,LLM)
!rt
!rt
!rt
      I00   = 0
      J00   = 0
!---------LREWND = .FALSE. IF A FULL YEAR OF HISTORY DATA IS NOT AVAILABLE!
      IF (NTEND-NTBEG .NE. 8760)   LREWND = .FALSE.
      WRITE(lun13,982) LHEADR,LREWND,NTBEG,NTEND,NAVG,NREAD
!xx
      II0 = 0
      JJ0 = 0
      LL0 = 0
!xx
      WRITE(lun13,983) II0,IIM,JJ0,JJM,LL0,LLM
!-------------CHECK WINDS & GRID SIZE-----------------------------------
      IF (IIM.LT.IM .OR. JJM.LT.JM) GOTO 91
      IF (IIM.NE.IMWIND .OR. JJM.NE.JMWIND .OR. LLM.NE.LMWIND) GOTO 91
!---------wind fields: TAUWB+NREAD = TAUW, winds = avg over NREAD hours,
!---------             pressure at end of intervval = TAUW
!---------locate wind field such that  TAUW .LE. ITAU and ITAU .LT. TAUW + NREAD
!---------map the location of ITAU into the range NTBEG to NTEND:-------
      ITAU1 = (ITAU /NREAD) *NREAD
!rt
!rt
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,906) number
        Write (lun13,900) ITAU,NREAD,ITAU1
      endif
!rt
!rt
!rt
      ITAU1 = MOD(ITAU1, 8760)
!rt
!rt
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,901) ITAU1
      endif
!rt
!rt
!rt
      IF (ITAU1 .LT. NTBEG)  ITAU1 = ITAU1 +((NTBEG -1) /8760 +1) *8760
!rt
!rt
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,902) ITAU1,NTBEG
      endif
!rt
!rt
!rt
      IF (ITAU1.GT.NTEND .OR. ITAU1.EQ.NTBEG) GOTO 91
      NSKIP = (ITAU1 -NTBEG) /NREAD - 1
!rt
!rt
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,903) NSKIP
      endif
!rt
!rt
!rt
!rt   DO ISKIP = 1, NSKIP
!rt     READ(3,ERR=93,END=95) TAUWB,TAUW
!rt   ENDDO
!rt   READ(3,ERR=93,END=95) TAUWB,TAUW,U,V,P,CWET,CDRY,BLH,EFFK,TAUX
!rt
!rt
!rt
      ihyper=NSKIP+1
      if (edit_type .eq. "debug") then
        Write (lun13,905) ihyper
      endif
!rt
      call met_readv (cdfid,ihyper,  &
     &   IWPAR,JWPAR,LWPAR,IIPAR,JJPAR,  &
     &   TAUWB,TAUW,TAUX,U,V,P,CWET,CDRY,BLH,EFFK)
!rt
      IF (TAUW .NE. TAUX)  GOTO 93
      NTAUW = (TAUW +0.01)
      ITAUW = MOD(NTAUW,8760)
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,904) ITAU,ITAUW,NTAUW,TAUW
      endif
!rt
!rt
!rt
      GOTO 99
!---------read in next wind/press/etc record----------------------------
!rt41 READ(3,ERR=93,END=61) TAUWB,TAUW,U,V,P,CWET,CDRY,BLH,EFFK,TAUX
!rt
!rt
!rt
 41   continue
      ihyper=ihyper+1
      if (ihyper.gt.number) goto 61
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,905) ihyper
      endif
!rt
      call met_readv (cdfid,ihyper,  &
     &   IWPAR,JWPAR,LWPAR,IIPAR,JJPAR,  &
     &   TAUWB,TAUW,TAUX,U,V,P,CWET,CDRY,BLH,EFFK)
!rt
      IF (TAUW .NE. TAUX)  GOTO 93
      NTAUW = (TAUW+0.1)
      ITAUW = MOD(NTAUW-1,8760) + 1
!rt
      if (edit_type .eq. "debug") then
        Write (lun13,904) ITAU,ITAUW,NTAUW,TAUW
      endif
!rt
!rt
!rt
      GOTO 99
!---------end of wind file, REWIND & recycle----------------------------
!rt61 IF (LREWND)   THEN
!rt     REWIND 3
!rt     READ(3,ERR=93,END=95) NTAU1,NTAU2
!rt     WRITE (lun13,'(A,2I10)') ' REWOUND GCM WIND FILE',NTAU1,NTAU2
!rt     READ(3,ERR=93,END=95) TAUWB,TAUW,U,V,P,CWET,CDRY,BLH,EFFK,TAUX
!rt
 61   if (LREWND) then
!rt
!rt
!rt
      call met_readh (cdfid,number,  &
     &   NTBEG,NTEND,NAVG,NREAD,  &
     &   I00,IIM,J00,JJM,L00,LLM)
!rt
      WRITE (lun13,'(A,2I10)') ' REWOUND GCM WIND FILE',NTBEG,NTEND
!rt
      ihyper=1
!rt
      call met_readv (cdfid,ihyper,  &
     &   IWPAR,JWPAR,LWPAR,IIPAR,JJPAR,  &
     &   TAUWB,TAUW,TAUX,U,V,P,CWET,CDRY,BLH,EFFK)
!rt
!rt
!rt
        IF (TAUW .NE. TAUX)  GOTO 93
        NTAUW = (TAUW+0.1)
        ITAUW = MOD(NTAUW-1,8760) + 1
       ELSE
         CALL EXITC('=====NOT ALLOWED TO REIWND WIND FILE====')
      ENDIF
      GOTO 99
!----------------------------------------------------------------------------
!---------error on initial wind files------------------------------------
   91 CALL EXITC('****ERROR ILLOGICAL WIND TAPE/TRACER RUN')
!---------error on unformatted input------------------------------------
   93 CALL EXITC(' =====ERROR ON WIND/CTM FILES=====      ')
!---------unexpected end of wind file-----------------------------------
   95 CALL EXITC(' =====UNEXPECTED EOF ON WINDS=====      ')
  982 FORMAT(' GCM DATA: LHEADR/LREWND/NTBEG/NTEND/NAVG/NREAD',2L2,4I10)
  983 FORMAT(' GCM DATA: II0/IIM/JJ0/JJM/LL0/LLM',6I5)
!rt
!rt
!rt
 900  format (1x,'ITAU=',i8,1x,'NREAD=',i8,1x,'ITAU1=',i8)
 901  format (1x,'ITAU1=',i8)
 902  format (1x,'ITAU1=',i8,1x,'NTBEG=',i8)
 903  format (1x,'NSKIP=',i8)
 904  format (1x,'WIND',18x,'ITAU=',i8,1x,'ITAUW=',i8,  &
     &   1x,'NTAUW=',i8,1x,'TAUW=',e13.6)
 905  format (1x,'WIND',1x,'ihyper=',i8)
 906  format (1x,'number=',i8)
!rt
!rt
!rt
   99 RETURN
      END



      SUBROUTINE EXITC (MESSAG) 4
!xx   IMPLICIT NONE
!xx
#     include "uci.h"
!xx
      CHARACTER*40 MESSAG
!xx$$
      save
!xx$$
      WRITE(lun13,'(A)') MESSAG
      STOP
      END



      subroutine jshift1(inarr,outarr) 1
!       shift mass/ad array in j by one index

#include "gmi_dims.h"
#include "uci.h"

      real*8 outarr(i1:i2, ju1:j2, k1:k2)
      real*8 inarr(IIPAR,JJPAR,LLPAR)
      integer  il, ij, ik

        do 888 ik=k1,k2
          do 887 ij=ju1,j2
            do 886 il=i1,i2
               outarr(il,ij,ik)=inarr(il,ij+1,ik)
  886       continue
  887     continue
  888   continue

      return
      end



      subroutine ad_to_mass(inarr,outarr)
!       convert ad array to mass Gmimod form array

#include "gmi_dims.h"
#include "uci.h"

      real*8 outarr(i1:i2, ju1:j2, k1:k2)
      real*8 inarr(IIPAR,JJPAR,LLPAR)
      integer  il, ij, ik

        do 888 ik=k1,k2
          do 887 ij=1,JJPAR
            do 886 il=i1,i2
               outarr(il,ij,ik)=inarr(il,ij,ik)
  886       continue
  887     continue
  888   continue

      return
      end



      subroutine jshiftn(inarr,outarr)
!     shift stt to const and convert to mixing ratio for diag0_gmi

#include "gmi_dims.h"
#include "uci.h"

      real*8 outarr(i1:i2, ju1:j2, k1:k2, 1:num_species)
      real*8 inarr(IIPAR,JJPAR,LLPAR,NNPAR)
      integer  il, ij, ik, ic

      do 889 ic=1,num_species
        do 888 ik=k1,k2
          do 887 ij=ju1,j2
            do 886 il=i1,i2
               outarr(il,ij,ik,ic)=tcvv(ic)*inarr(il,ij+1,ik,ic)/  &
     &                                           ad(il,ij+1,ik)
  886       continue
  887     continue
  888   continue
  889 continue

      return
      end



      subroutine stt_to_const(inarr,outarr)
!     routine to convert stt to const form for diagnostics, no shift
!     for gmi_diag_save

#include "gmi_dims.h"
#include "uci.h"

      real*8 outarr(i1:i2, ju1:j2, k1:k2, 1:num_species)
      real*8 inarr(IIPAR,JJPAR,LLPAR,NNPAR)
      integer  il, ij, ik, ic

      do 889 ic=1,num_species
        do 888 ik=k1,k2
          do 887 ij=1,JJPAR  ! FIX ME
            do 886 il=i1,i2
               outarr(il,ij,ik,ic)=tcvv(ic)*inarr(il,ij,ik,ic)/  &
     &                                           ad(il,ij,ik)
  886       continue
  887     continue
  888   continue
  889 continue

      return
      end


!-----------------------------------------------------------------------------
!
! ROUTINE
!   Get_Diag_1d_Dump_Flag
!
! DESCRIPTION
!   This routine determines whether or not it is time for the UCI 1D
!   diagnostics to be dumped.
!
! ARGUMEMTS
!   dynmd     : year/month/day (YYYYMMDD) at end of current time step
!   dump_diag : time to dump diagnostics?
!
!-----------------------------------------------------------------------------

!c    subroutine Get_Diag_1d_Dump_Flag (dnymd, dump_diag)

!c    implicit none

!c#     include "gmi_control.h"
!c#     include "gmi_gen_io.h"


!     ----------------------
!     Argument declarations.
!     ----------------------

!c    integer :: dnymd
!c    logical :: dump_diag


!     ----------------------
!     Variable declarations.
!     ----------------------

!c    logical, save :: first_diag_uci = .true.

!c    integer :: idumday
!c    integer :: idumyear
!c    integer :: this_month_diag_uci

!c    integer, save :: save_month_diag_uci = 0


!     ----------------
!     Begin execution.
!     ----------------

!c    dump_diag = .false.

!     ===============
!c    if (pr_diag_1d) then
!     ===============

!c      if (pr_diag_uci_step_interval < 0) then

!         -------------------------------------------
!         UCI diagnostic output at monthly intervals.
!         -------------------------------------------

!c        if (first_diag_uci) then

!c          dump_diag      = .true.

!c          call Split_Time_Flds
!c   &        (dnymd, idumyear, save_month_diag_uci, idumday)

!c          first_diag_uci = .false.

!c        else

!c          call Split_Time_Flds
!c   &        (dnymd, idumyear, this_month_diag_uci, idumday)

!c          if (this_month_diag_uci /= save_month_diag_uci) then

!c            dump_diag           = .true.
!c            save_month_diag_uci = this_month_diag_uci

!c          end if

!c        end if

!       ----------------------------------------------------------
!       The +1 below is because the clock has not been ticked yet.
!       ----------------------------------------------------------

!c      else if (Mod (num_time_steps+1,
!c   &                pr_diag_uci_step_interval) == 0) then

!         --------------------------------------------
!         UCI diagnostic output at specified interval.
!         --------------------------------------------

!c        dump_diag = .true.

!c      end if

!c    end if


!c    return

!c    end