MODULE output_module
   
   USE input_module
   USE module_model_basics
   USE module_arrays
   
   integer                             :: time, rec, ivars   !!! process time ; rec in .dat file ; variables in .dat file
   integer                             :: cunit, dunit       !!! .ctl and .dat output files
   character (len=128)                 :: ctlfile, datfile, v5dfile
   character (len=200), dimension(200) :: ctl_var_string     !!! List of fields foe .ctl file
   character (len=2000)                :: could_not_find     !!! Keep list of what we wrote out
   character (len=19)                  :: tdef_date          !!! Output start date
   character (len=40)                  :: tdef               !!! tdef string in .ctl file
   character (len=10), dimension(100)  :: varnames 
   integer                             :: numvars

   CONTAINS

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: create_ctl     
   ! Purpose: Create the .ctl file                  
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE create_ctl ()

      IMPLICIT NONE
  
      ! Local variables
      real, dimension(4)               :: xlat_a, xlon_a
      real                             :: dxll, slack, temp
      real                             :: abslatmin, abslonmax, abslonmin, abslatmax
      integer                          :: i, k, ipoints, jpoints, ilon 
      integer                          :: ipole
      real                             :: xi, r, re, xj, dx_ps, earthr, radpd
      real                             :: rlat, rlong, ref_lon, ref_lat
      integer                          :: file_start, file_end
      integer                          :: we_dim, sn_dim


      
      file_start = INDEX(datfile,"/",BACK=.TRUE.)
      file_end   = len_trim(datfile)
      write(cunit,'("dset ^",A)') datfile(file_start+1:file_end)
#ifdef bytesw
   write (cunit, '("options byteswapped")')
#endif
      write(cunit, '("undef 1.e37")')

      IF ( output_title /= '   ') THEN
        write(cunit,'("title ",A)') trim(output_title)
      ELSE
        write(cunit,'("title ",A)') trim(title)
      END IF

      IF (map_proj == 1) THEN                                     !! Lambert Projection


!!!!!!!
!!!!!!! TEMP
!!!!!!!
IF ( mercator_defs ) THEN
   WRITE(cunit,'("xdef ",i4," levels ")') west_east_dim
     DO i = 1,west_east_dim
        WRITE(cunit,*) XLONG(i,int(west_east_dim/2))
     ENDDO
   WRITE(cunit,'("ydef ",i4," levels ")') south_north_dim
     DO i = 1,south_north_dim
        WRITE(cunit,*) XLAT(int(south_north_dim/2),i)
     END DO
ELSE

        ! Make sure truelat1 is the larger number
         IF (truelat1 < truelat2) THEN
            temp = truelat1
            truelat1 = truelat2
            truelat2 = temp
         END IF

         !!WRITE(cunit,'("pdef ",i4," ",i3," lcc ",f9.5," ",f10.5," ", &   !! rounding off
         WRITE(cunit,'("pdef ",i4," ",i3," lcc ",f7.3," ",f8.3," ", &    
&             f8.3," ",f8.3," ",f9.5," ",f9.5," ",f10.5," ",         &
&             f10.3," ",f10.3)')                                       &
              west_east_dim,south_north_dim,cen_lat,cen_lon,         &
              (west_east_dim+1)/2.,(south_north_dim+1)/2.,           &
              truelat1,truelat2,stand_lon,dx,dy


         xlon_a(1) = XLONG(1,1)
         xlon_a(2) = XLONG(1,south_north_dim)
         xlon_a(3) = XLONG(west_east_dim,1)
         xlon_a(4) = XLONG(west_east_dim,south_north_dim)

         !check for dateline
         ilon = 0
         if ( abs(xlon_a(1) - xlon_a(3)) .GT. 180. ) ilon = 1
         if ( abs(xlon_a(2) - xlon_a(4)) .GT. 180. ) ilon = 1

         abslatmin = minval(XLAT)
         abslatmax = maxval(XLAT)
         abslonmin =  99999.
         abslonmax = -99999.
         DO i=1,4
           IF ( xlon_a(i) .lt. 0.0 .AND. ilon .eq. 1 ) THEN
             abslonmin=min(abslonmin,360.+xlon_a(i))
             abslonmax=max(abslonmax,360.+xlon_a(i))
           ELSE
             abslonmin=min(abslonmin,xlon_a(i))
             abslonmax=max(abslonmax,xlon_a(i))
           ENDIF
         END DO

!         dxll = (dx/1000.)/111./2.
         dxll = (dx/1000.)/59./2.
         slack = ((dx/1000.)*3.)/100.
         ipoints = int(((abslonmax-abslonmin)+(3*slack))/dxll)
         jpoints = int(((abslatmax-abslatmin)+(3*slack))/dxll)

         WRITE(cunit,'("xdef ",i4," linear ",f10.5," ",f12.8)') ipoints, &
                      (abslonmin-1.5*slack),dxll
         WRITE(cunit,'("ydef ",i4," linear ",f10.5," ",f12.8)') jpoints, &
                      (abslatmin-1.5*slack),dxll

!!!!!!!!!!!!!
!!!!!!!!!!!!! TEMP
!!!!!!!!!!!!!
ENDIF


      ENDIF

      IF (map_proj == 2) THEN                                     !! Polar Stereo Projection

         ipole = 1
         IF (truelat1 .lt. 0.0) ipole = -1
         radpd  = .01745329
!         earthr = 6.3712E6
	 earthr = 3.3972E6

         dx_ps = dx/( (1.0+sin(ipole*TRUELAT1*radpd))/(1.0+sin(60*radpd)) ) !! true dx (meter) at 60degrees
         re    = ( earthr * (1.0 + sin(60*radpd)) ) / dx_ps          !! 


         IF ( ipole == 1 ) THEN
           rlong  = (XLONG(1,1) - stand_lon) * radpd
           rlat   = XLAT(1,1) * radpd

           r    =  (re * cos(rlat)) / (1.0 + sin(rlat))
           xi   =  (1. - r * sin(rlong))
           xj   =  (1. + r * cos(rlong))
           
           write(cunit,'("pdef ",i3," ",i3," nps ",f8.3," ",f8.3," ", &
&                f12.4," ",f12.7)') west_east_dim, south_north_dim,   &
                 xi, xj, stand_lon, dx_ps*0.001

         ELSE IF ( ipole == -1 ) THEN
           rlong  = (180.0 + XLONG(1,1) - stand_lon) * radpd
           rlat   = ipole*XLAT(1,1) * radpd

           r    =  (re * cos(rlat)) / (1.0 + sin(rlat))
           xi   =  (1. - ipole * r * sin(rlong))
           xj   =  (1. + r * cos(rlong))

           write(cunit,'("pdef ",i3," ",i3," sps ",f8.3," ",f8.3," ", &
&                f12.4," ",f12.7)') west_east_dim, south_north_dim,   &
                 xi, xj, (180+stand_lon), -0.001*dx_ps
         END IF


         abslonmin = minval(XLONG)
         abslonmax = maxval(XLONG)
         abslatmin = minval(XLAT)
         abslatmax = maxval(XLAT)
   
!!         dxll = (dx_ps/1000.)/111./2.
         dxll = (dx_ps/1000.)/59./2.
         ipoints = int(((abslonmax-abslonmin)+2.)/dxll)-2
         jpoints = int(((abslatmax-abslatmin)+2.)/dxll)-3

         ref_lon = abslonmin-1.
         IF ( abslonmax .ge. 175. .AND. abslonmax .le. 180.) THEN
           IF ( abslonmin .ge. -180. .AND. abslonmin .le. -175.) THEN
             ref_lon = 0.0
           END IF
         END IF

         ref_lat = abslatmin-1.
         !!IF ( ipole == -1 ) ref_lat = abslatmax+1.
         IF ( ipole == -1 ) ref_lat = abslatmin
   
         WRITE(cunit,'("xdef ",i4," linear ",f6.1," ",f12.8)') ipoints, &
                      ref_lon,dxll
         WRITE(cunit,'("ydef ",i4," linear ",f6.1," ",f12.8)') jpoints, &
                      ref_lat,dxll

      ENDIF

      IF (map_proj == 3) THEN                         !! Mercator Projection
   
          !check for dateline
          ilon = 0
          IF ( abs(XLONG(1,1) - XLONG(west_east_dim,1)) .GT. 180. ) ilon = 1
          IF ( abs(XLONG(1,south_north_dim) - XLONG(west_east_dim,south_north_dim)) .GT. 180. ) ilon = 1
   
          IF ( mercator_defs .AND. map_proj == 3 ) THEN
            WRITE(cunit,'("xdef ",i4," levels ")') west_east_dim
            DO i = 1,west_east_dim
              WRITE(cunit,*) XLONG(i,1)
              !!IF ( XLONG(i,1) .lt. 0.0 ) THEN
                !!WRITE(cunit,*) XLONG(i,1) + 180.0
              !!ELSE
                !!WRITE(cunit,*) XLONG(i,1)
              !!END IF
            END DO
          ELSE
            IF ( ilon == 1 ) THEN
              WRITE(cunit,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
                    west_east_dim,XLONG(1,1), &
                    (abs(XLONG(1,1)-(360.+XLONG(west_east_dim,south_north_dim)))/(west_east_dim-1))
            ELSE
              WRITE(cunit,'("xdef ",i4," linear ",f9.4," ",f8.4)')   &
                    west_east_dim,XLONG(1,1), &
                    (abs(XLONG(1,1)-XLONG(west_east_dim,south_north_dim))/(west_east_dim-1))
            ENDIF
          ENDIF

          IF ( mercator_defs .AND. map_proj == 3 ) THEN
             WRITE(cunit,'("ydef ",i4," levels ")') south_north_dim
             DO i = 1,south_north_dim
               WRITE(cunit,*) XLAT(1,i)
             END DO
          ELSE
             WRITE(cunit,'("ydef ",i4," linear ",f9.4," ",f8.4)')   &
                       south_north_dim,XLAT(1,1),(abs(XLAT(1,1)-XLAT(west_east_dim,south_north_dim))/(south_north_dim-1))
          ENDIF

      ENDIF


      IF (map_proj == 0) THEN                         !! Ideal
        we_dim = west_east_dim
        sn_dim = south_north_dim
        IF (we_dim == 2) we_dim = 1
        IF (sn_dim == 2) sn_dim = 1
        write (cunit, '("xdef  ",i4, " linear 0 0.0001")') we_dim
        write (cunit, '("ydef  ",i4, " linear 0 0.0001")') sn_dim
      ENDIF
      !! END of projection specific section

      IF (iprogram .ge. 6 .and. (vertical_type=='p'.or.vertical_type=='z')) THEN
        WRITE(cunit,'("zdef  ",i3, " levels  ")') number_of_zlevs
        DO k = 1,number_of_zlevs
           WRITE(cunit,'(f10.5)') interp_levels(k)
        END DO
      ELSE
        WRITE(cunit,'("zdef ",i4," linear 1 1  ")') bottom_top_dim
      ENDIF

       IF (tdef_date == '0000-00-00_00:00:00') THEN
         tdef = '   1 linear 00z01jan2000  1hr'
       ELSE
         IF ( output_type == 'grads' ) THEN
                CALL tdef_calc_grads () 
         ELSE 
                CALL tdef_calc ()
         ENDIF
       END IF
       WRITE(cunit,'("tdef ",A)') tdef

       WRITE(cunit,'("VARS  ",i3)') ivars
       DO i = 1,ivars
          WRITE(cunit,'(A)') trim(ctl_var_string(i))
       ENDDO
       WRITE(cunit,'("ENDVARS")')
       
       DO i=1,iatts
         WRITE(cunit,'(A)') trim(catts(i))
       ENDDO


   END SUBROUTINE create_ctl

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: tdef_calc      
   ! Purpose: Make the tdef line for the .ctl file
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE tdef_calc ()

     IMPLICIT NONE

     character (len=3)              :: cmonth
     integer                        :: year,month,day,hour
  
     tdef = '                               '
  
     !! Time comes in as YYYY-MM-DD_HH:MM:SS
     !!                  1234 67 90 23 45 89
     READ(tdef_date(1:4),*)   year
     READ(tdef_date(6:7),*)   month
     READ(tdef_date(9:10),*)  day
     READ(tdef_date(12:13),*) hour

    IF (month == 1) cmonth = 'M1'
    IF (month == 2) cmonth = 'M2'
    IF (month == 3) cmonth = 'M3'
    IF (month == 4) cmonth = 'M4'
    IF (month == 5) cmonth = 'M5'
    IF (month == 6) cmonth = 'M6'
    IF (month == 7) cmonth = 'M7'
    IF (month == 8) cmonth = 'M8'
    IF (month == 9) cmonth = 'M9'
    IF (month ==10) cmonth = 'M10'
    IF (month ==11) cmonth = 'M11'
    IF (month ==12) cmonth = 'M12'

    IF ( year < 100 ) year = 2000

    WRITE (tdef,'(i4," linear ",i2.2,"Z",i2.2,A,i4,"  ",i6,"MN")') time, hour, day, cmonth, year, interval_seconds

   END SUBROUTINE tdef_calc

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: tdef_calc_grads      
   ! Purpose: Make the tdef line for the .ctl file
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE tdef_calc_grads ()

     IMPLICIT NONE

     character (len=3)              :: cmonth
     integer                        :: year,month,day,hour

     tdef = '                               '

     !! Time comes in as YYYY-MM-DD_HH:MM:SS
     !!                  1234 67 90 23 45 89
     READ(tdef_date(1:4),*)   year
     READ(tdef_date(6:7),*)   month
     READ(tdef_date(9:10),*)  day
     READ(tdef_date(12:13),*) hour

    IF (month == 1) cmonth = 'jan'
    IF (month == 2) cmonth = 'fev'
    IF (month == 3) cmonth = 'mar'
    IF (month == 4) cmonth = 'avr'
    IF (month == 5) cmonth = 'may'
    IF (month == 6) cmonth = 'jun'
    IF (month == 7) cmonth = 'jul'
    IF (month == 8) cmonth = 'aug'
    IF (month == 9) cmonth = 'sep'
    IF (month ==10) cmonth = 'oct'
    IF (month ==11) cmonth = 'nov'
    IF (month ==12) cmonth = 'dec'

    IF ( year < 100 ) year = 2000

    WRITE (tdef,'(i4," linear ",i2.2,"Z",i2.2,A,i4,"  ",i6,"MN")') time, hour, day, cmonth, year, interval_seconds

   END SUBROUTINE tdef_calc_grads


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: write_dat
   ! Purpose: Write the data to the .dat file. Also keep record of what has been
   !          written for creation of the .ctl file
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   SUBROUTINE write_dat (data_out, nx, ny, nz, cname, cdesc, cunits)

      IMPLICIT NONE

      ! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS
      include 'v5df.h'

      ! Arguments
      integer                        :: nx, ny, nz
      real, dimension(nx,ny,nz)      :: data_out
      character (len=128)            :: cname, cunits, cdesc

      ! Local variables
      integer                        :: ii, jj, kk
      integer                        :: is_there, ret, ivar
      character (len=20)             :: dummy
      real, dimension(ny,nx,nz)      :: v5d_data_out


!     IF ( output_type == 'ncdf' ) THEN   !!! Create the netcdf file
!        IF ( nx == 2 ) nx = 1
!        IF ( ny == 2 ) ny = 1
!        !DO kk=1,nz
!        !  rec = rec + 1
!        !  WRITE(dunit,rec=rec) ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny)
!        !END DO
!        CALL writencdf(unite nx,ny,nz,"u","W.m-2",3,data_out)
!     END IF



      IF ( output_type == 'grads' ) THEN   !!! Create the grads file
        IF ( nx == 2 ) nx = 1
        IF ( ny == 2 ) ny = 1
        DO kk=1,nz
          rec = rec + 1
          WRITE(dunit,rec=rec) ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny)
        END DO
      END IF

!!!!****MARS IDL
      IF ( output_type == 'idl' ) THEN   !!! Create the IDL ascii file
        IF ( nx == 2 ) nx = 1
        IF ( ny == 2 ) ny = 1
        DO kk=1,nz 
          rec = rec + 1
!!!****MARS
!          WRITE(dunit,rec=rec) ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny)
          WRITE(dunit,FMT='(F20.8)') ((data_out(ii,jj,kk),ii=1,nx),jj=1,ny)
        END DO
      END IF



#ifdef V5D
      IF ( output_type == 'v5d' ) THEN   !!! Create the vis5d file
        DO kk = 1,nz
        DO jj = 1,nx
        DO ii = 1,ny
          v5d_data_out(ii,jj,kk) = data_out(jj,ny-ii+1,kk)
        END DO
        END DO
        END DO


        ivar = 0
        DO ii = 1,numvars
          IF ( trim(varnames(ii)) == trim(cname) ) THEN 
            ivar = ii
          END IF
        END DO 
        ret = v5dwrite ( time+1, ivar, v5d_data_out )
      END IF
#endif


      IF (time == 0) THEN
        IF ( cname(1:3) == 'clf' ) then
          dummy = ",clfr,"
        ELSE
          dummy = ","//trim(cname)//","
        END IF
        is_there = INDEX(could_not_find,trim(dummy))
        DO WHILE ( is_there /= 0 )
          could_not_find = could_not_find(1:is_there)//could_not_find(is_there+len_trim(dummy):len_trim(could_not_find))
          is_there = INDEX(could_not_find,trim(dummy))
        END DO
        ivars = ivars + 1
        WRITE(ctl_var_string(ivars),'(A,"  ",i4,"  0  ",A," (",A,")")') cname(1:10),nz,trim(cdesc),trim(cunits)
      ENDIF


   END SUBROUTINE write_dat

!   SUBROUTINE writencdf(unite,nx,ny,nz,"u","W.m-2",3,data_out)
!
!
!        ! NETCDF file
!        fichnom="diagfi.nc"
!        if (firstnom.eq.'1234567890') then
!        firstnom = nom
!        ierr = NF_CREATE  (fichnom, NF_CLOBBER, nid)
!        ! Build dimensions
!        ierr = NF_DEF_DIM (nid, "Time",      NF_UNLIMITED, idim)
!        ierr = NF_DEF_DIM (nid, "latitude",  ny,           idim_lat)
!        ierr = NF_DEF_DIM (nid, "longitude", nx,           idim_lon)
!        ierr = NF_DEF_DIM (nid, "altitude",  nz,           idim_alt)
!        ! Write dimensions attributes
!        ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim, varid)
!        ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",  4, "Time")
!        ierr = NF_PUT_ATT_TEXT (nid, varid, 'units',     29, "Time")
!        ierr = NF_ENDDEF(nid)
!        !
!        ierr = NF_DEF_VAR (nid, "latitude",  NF_FLOAT,  1, idim_lat, varid)
!        ierr = NF_PUT_ATT_TEXT (nid, varid, 'units',     13, "degrees_north")
!        ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 14, "North latitude")
!        ierr = NF_ENDDEF(nid)
!        ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT,  1, idim_lon, varid)
!        ierr = NF_PUT_ATT_TEXT (nid, varid, 'units',     12, "degrees_east")
!        ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 14, "East longitude")
!        ierr = NF_ENDDEF(nid)
!        ierr = NF_DEF_VAR (nid, "altitude",  NF_FLOAT,  1, idim_alt, varid)
!        ierr = NF_PUT_ATT_TEXT (nid, varid, 'units',      2, "km")
!        ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name", 10, "pseudo-alt")
!        ierr = NF_PUT_ATT_TEXT (nid, varid, 'positive',   2, "up")
!        ierr = NF_ENDDEF(nid)
!        else
!        ierr = NF_OPEN(fichnom,NF_WRITE,nid)
!        endif
!        ! Advance time
!        ierr = NF_INQ_VARID (nid, "Time", varid)
!        ierr = NF_PUT_VARA_REAL (nid, varid, ntime, 1, date)
!        ! Write field
!        ierr = NF_INQ_VARID(nid,nom,varid)
!        if (ierr /= NF_NOERR) then
!        ! Create variable if does not exist
!        ierr = NF_INQ_DIMID(nid,"longitude",id(1))
!        ierr = NF_INQ_DIMID(nid,"latitude",id(2))
!        ierr = NF_INQ_DIMID(nid,"altitude",id(3))
!        ierr = NF_INQ_DIMID(nid,"Time",id(4))
!        write (*,*) "====================================="
!        write (*,*) "NETCDF: create ",nom
!        call def_var (nid,nom,titre,unite,4,id,varid,ierr)
!        endif
!   END SUBROUTINE writencdf


!!==================================================================================
!!==================================================================================
!
!subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)
!
!implicit none
!
!include "netcdf.inc"
!
!character (len=*) :: title,units,name
!integer :: nid,nbdim,nvarid,ierr
!integer, dimension(nbdim) :: dim
!
!ierr=NF_REDEF(nid)
!#ifdef NC_DOUBLE
!ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
!#else
!ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
!#endif
!if(ierr/=NF_NOERR) then
!   write(*,*) NF_STRERROR(ierr)
!   stop "in def_var"
!endif
!ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title",
!len_trim(adjustl(title)),adjustl(title))
!if(ierr/=NF_NOERR) then
!   write(*,*) NF_STRERROR(ierr)
!   stop "in def_var"
!endif
!ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units",
!len_trim(adjustl(units)),adjustl(units))
!if(ierr/=NF_NOERR) then
!   write(*,*) NF_STRERROR(ierr)
!   stop "in def_var"
!endif
!ierr = NF_ENDDEF(nid)
!
!end subroutine def_var


END MODULE output_module
