!WRF:PACKAGE:IO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      MODULE MAP_GOCART_DMS_V01

!
! Version 0.1 of the emissions standard initialization routine.
!
! This program takes formatted output of em99v3 NEI99v3 inventory and
! grids to a different grid.  Simple grid dumping done here.
!
!  Steven Peckham 1/25/08
!
! compile with:
!
! ifort -FR -convert big_endian map_gocart_dms_d01.F -I /opt/netcdf/3.6.1_intel_9.1/include -L/opt/netcdf/3.6.1_intel_9.1/lib -lnetcdf
!
! pgf90 -w -byteswapio -Mfree -Mlfs -I/opt/netcdf-3.5.0-pgi-5.2-4/include -L/opt/netcdf-3.5.0-pgi-5.2-4/lib -lnetcdf map_gocart_dms_d01.F
!
!-----------------------------------------------------------------------
! 
!  Fields to set before running:
! 
!  zfa                    elevation at grid cell top (m)
!  ix2                    x-dimension of output  data 
!  jx2                    y-dimension of output  data 
!  kx                     z-dimension of output  data 
!  DX                     horizontal grid spacing (m)
!  REBY                   Earth radius (km)
!  LAT1                   SW latitude (1,1) in degrees (-90->90N)
!  LON1                   SW longitude (1,1) in degrees (-180->180E)
!  STDLON                 standard longitude parallel to y-axis (-180->180E)
!  TRUELAT1               Northern most reference latitude of projection
!  TRUELAT2               Southern most reference latitude of projection (TRUELAT1 > TRUELAT2)
!  KNOWNI                 Origin point of domain, x-location
!  KNOWNJ                 Origin point of domain, y-location
!  HEMI                   1 for Northern Hemisphere, -1 for Southern Hemisphere
! 
!  DATADIR                  directory of emissions input data
!
!-----------------------------------------------------------------------

    use netcdf
     implicit none
!
!     Constants
      REAL       :: rad_per_deg = 0.0174533
      REAL       :: deg_per_rad = 57.2958

!
! - Input data dimensions for the GLOBAL DATA
!
      integer       ::  lon, ilon, lat, ilat, klev, iklev
      PARAMETER( lon=288, lat=181, klev=55)
      PARAMETER(ilon=289,ilat=182,iklev=56)

      integer                        :: status, system


! - Output data dimensions for the WRF domain(user specified)
!
      integer       ::  ix2,jx2,kx,kp
!     PARAMETER(IX2=39,JX2=39,KX=19,KP=KX+1)
      PARAMETER(IX2=139,JX2= 94,KX=40,KP=KX+1)
!
      character (len=80)   ::  dir_in='/p20/ap-fc/rt-aq/Chile/GOCART_DMS'
      character (len=80)   ::  fname_in_1='dms_1x1.25.nc'
      character (len=80)   ::  fname_out='wrfchem_gocartdms_raw.nc'

!
! Map information for output grid 
!          REBY = Earth radius (km)
!          DX = horizontal grid spacing (m)
!          LAT1 = SW latitude (1,1) in degrees (-90->90N)
!          LON1 = SW longitude (1,1) in degrees (-180->180E)
!          STDLON = standard longitude parallel to y-axis (-180->180E)
!          TRUELAT1 = Northern most reference latitude of projection
!          TRUELAT2 = Southern most reference latitude of projection (TRUELAT1 > TRUELAT2)
!          KNOWNI = origin point of domain, x-location
!          KNOWNJ = origin point of domain, y-location
!          HEMI   = 1 for Northern Hemisphere, -1 for Southern Hemisphere
!
      real       :: reby = 6370.
      real       :: dx = 36.E3
      real       :: lat1 =  38.00
      real       :: lon1 = -097.00
      real       :: stdlon = -97.0
      real       :: knowni = 70.0
      real       :: knownj = 47.0
      real       :: truelat1 = 38.001
      real       :: truelat2 = 38.000
      integer    :: hemi = +1

      CONTAINS

      SUBROUTINE  GET_GOCART_DMS_DATA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      include 'netcdf.inc'

      real, dimension(lat)                     :: lat_g
      real, dimension(lon)                     :: lon_g
      real, dimension(12)                      :: time_g

      real, dimension(lon, lat, 12)            :: dms_g
      real, dimension(lon, lat)            :: dms_m
      real, dimension(ix2, jx2, 12)            :: dms0

      real, dimension(ix2,jx2)             :: xlat,xlon, dms_out
      real                             :: dist, a1, a2, b1, b2

      integer                        :: lat_id,lon_id
      integer                        :: time_id, dms_id

      integer     ::  iloc, jloc
      integer     ::  i, j, k, jmx
      integer     ::  i2, j2
      real        ::  x2, y2, fac
      real        ::  top, bot
      integer     ::  kw, kg
      integer     ::  kbot, ktop, ! month

  ! This will be the netCDF ID for the file and data variable.
  integer :: ncid, varid

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      print *,' OPEN: '
! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file.
!     call check( nf_open(TRIM(dir_in)//TRIM(fname_in_1), NF_NOWRITE, ncid) )
      call check( nf_open(TRIM(fname_in_1), NF_NOWRITE, ncid) )

      print *,' LAT: '
! Get the varid of the data variable, based on its name.
      call check( nf_inq_varid(ncid, "lat", lat_id) )

! Read the data.
      call check( nf_get_var_real(ncid, lat_id, lat_g) )

      print *,' LON: '
! Get the varid of the data variable, based on its name.
      call check( nf_inq_varid(ncid, "lon", lon_id) )

! Read the data.
      call check( nf_get_var_real(ncid, lon_id, lon_g) )

      print *,' TIME: '
! Get the varid of the data variable, based on its name.
      call check( nf_inq_varid(ncid, "time", time_id) )

! Read the data.
      call check( nf_get_var_real(ncid, time_id, time_g) )

      print *,' DMS: ',time_g(1:12)
! Get the varid of the data variable, based on its name.
      call check( nf_inq_varid(ncid, "DMSO", dms_id) )

! Read the data.
      call check( nf_get_var_real(ncid, dms_id, dms_g) )

! Close the file, freeing all resources.
      call check( nf_close(ncid) )

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!     Fill xlat and xlon arrays for the WRF domain
      do j=1,jx2
      do i=1,ix2
           CALL MAPCF(float(i),float(j),xlat(i,j),xlon(i,j))
      enddo
      enddo
!     print *,' LAT_WRF: ', xlat(1,:)
!     print *,' LON_WRF: ', xlon(:,1)

!     Find where xlat and xlon arrays for the WRF domain are located on global domain
!     do month=1,12
!     dms_m(:,:)=dms_g(:,:,month)
!     write(20+month)dms_m 
!     enddo
      do j=1,lat-1
      do i=1,lon-1
!          print *,i,j,lat_g(j),lon_g(i)

! wrf domain loop
      iloc = 0
      jloc = 0
      do j2=1,jx2   !jx2
      do i2=1,ix2   !ix2
!          Find southwestern corner (lower left side of glocal grid box)
           if(  ( xlat(i2, j2) >= lat_g(j) .AND. xlat(i2,j2) <= lat_g(j+1)  ) .AND.  &
                ( xlon(i2, j2) >= lon_g(i) .AND. xlon(i2,j2) <= lon_g(i+1)  ) ) THEN
              iloc = i
              jloc = j

!          print *,' LAT_G: ', lat_g(jloc),lat_g(jloc+1)
!          print *,' LON_G: ', lon_g(iloc),lon_g(iloc+1)
!          print *,' WRF: ', xlat(i2,j2), xlon(i2,j2)

           a1 = ( xlat(i2,j2) - lat_g(jloc)) / (lat_g(jloc+1) - lat_g(jloc))
           a2 = 1. - a1
           b1 = ( xlon(i2,j2) - lon_g(iloc)) / (lon_g(iloc+1) - lon_g(iloc))
           b2 = 1. - b1

           dms_out(i2,j2) = b2*a2*dms_g(iloc  ,jloc  ,7) + b1*a2*dms_g(iloc+1,jloc  ,7) &
                          + b2*a1*dms_g(iloc  ,jloc+1,7) + b1*a1*dms_g(iloc+1,jloc+1,7)
 
!          return
!          print *, dms_out(i2,:,j2)
           endif   ! if test to find SW corner

      enddo  ! end wrf domain loops
      enddo

      enddo ! end global domain
      enddo


! Open output files
      OPEN(12,FILE='wrf_gocart_dms.txt',FORM='FORMATTED')
      write(12,'(12E15.5)') dms_out(:,:)

      CLOSE(12)

      END SUBROUTINE  GET_GOCART_DMS_DATA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      SUBROUTINE MAPCF (XI,YJ,XLATP,XLONP)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!                                                                     C
!                                                                     C
!     THIS SUBROUTINE COMPUTES THE LATITUDE AND LONGITUDE FROM MODEL  C
!     INDEXES OF A POINT.                                             C
!                                                                     C
!     INPUT :                                                         C
!                                                                     C
!        XI : X COORDINATE OF THE POINT IN MODEL INDEX.               C
!                                                                     C
!        YJ : Y COORDINATE OF THE POINT IN MODEL INDEX.               C
!                                                                     C
!     OUTPUT :                                                        C
!                                                                     C
!        XLATP : LATITUDE OF THE POINT.                               C
!                                                                     C
!        XLONP : LONGITUDE OF THE POINT.                              C
!                                                                     C
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
      REAL :: XI, YJ, XLATP, XLONP
      REAL :: rebydx
      REAL :: chi1, chi2, chi
      INTEGER :: jnew, inew, i, j
      REAL :: reflon, ala1, alo1, scale_top
      REAL :: deltalon1,cone, tl1r, ctl1r
      REAL :: rsw, arg, polei, polej, xx, yy, r2, r

      rebydx = reby/dx*1e3

      chi1 = (90. -hemi*truelat1)*rad_per_deg
      chi2 = (90. -hemi*truelat2)*rad_per_deg

      inew = hemi * int(xi)
      jnew = hemi * int(yj)

      reflon = stdlon + 90.
      ala1 = lat1 * rad_per_deg
      alo1 = (lon1 - reflon) * rad_per_deg
      scale_top = 1. + hemi * SIN(truelat1 * rad_per_deg)

      deltalon1 = lon1 - stdlon
      IF (deltalon1 .gt. +180.) deltalon1 = deltalon1 - 360.
      IF (deltalon1 .lt. -180.) deltalon1 = deltalon1 + 360.

      IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
         cone = ALOG10(COS(truelat1*rad_per_deg)) - &
                ALOG10(COS(truelat2*rad_per_deg))
         cone=cone/(ALOG10(TAN((45.0-ABS(truelat1)/2.0)*rad_per_deg))- &
                ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
      ELSE
         cone = SIN(ABS(truelat1)*rad_per_deg )
      ENDIF

      tl1r = truelat1 * rad_per_deg
      ctl1r = COS(tl1r)
      rsw = rebydx * ctl1r/cone * &
             (TAN((90.*hemi-lat1)*rad_per_deg/2.) / &
              TAN((90.*hemi-truelat1)*rad_per_deg/2.))**cone

      arg = cone*(deltalon1*rad_per_deg)
      polei = hemi*knowni - hemi * rsw * SIN(arg)
      polej = hemi*knownj + rsw * COS(arg)

      xx = inew - polei
      yy = polej - jnew
      r2 = (xx*xx + yy*yy)
      r = SQRT(r2)/rebydx

   ! Convert to lat/lon
      IF (r2 .EQ. 0.) THEN
         xlatp = hemi * 90.
         xlonp = stdlon
      ELSE
         ! Longitude
         xlonp = stdlon + deg_per_rad * ATAN2(hemi*xx,yy)/cone
         xlonp = AMOD(xlonp+360., 360.)

         ! Latitude.  Latitude determined by solving an equation adapted
         ! from:
         !  Maling, D.H., 1973: Coordinate Systems and Map Projections
         ! Equations #20 in Appendix I.

         IF (chi1 .EQ. chi2) THEN
            chi = 2.0*ATAN( ( r/TAN(chi1) )**(1./cone) * TAN(chi1*0.5) )
         ELSE
            chi = 2.0*ATAN( (r*cone/SIN(chi1))**(1./cone)*TAN(chi1*0.5))
         ENDIF
          xlatp = (90.0-chi*deg_per_rad)*hemi
      ENDIF

      IF (xlonp .GT. +180.) xlonp = xlonp - 360.
      IF (xlonp .LT. -180.) xlonp = xlonp + 360.

!     print *, xlonp,xlatp
      RETURN
!
      END SUBROUTINE  MAPCF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


      SUBROUTINE LAMBC (XLATP,XLONP,XI,YJ)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                     C
!                                                                     C
!     THIS SUBROUTINE COMPUTES THE MODEL INDEXES OF A POINT FROM      C
!     LATITUDE AND LONGITUDE.                                         C
!                                                                     C
!     INPUT :                                                         C
!                                                                     C
!        XLATP : LATITUDE OF THE POINT.                                C
!                                                                     C
!        XLONP : LONGITUDE OF THE POINT.                               C
!                                                                     C
!     OUTPUT :                                                        C
!                                                                     C
!        XI : X COORDINATE OF THE POINT IN MODEL INDEX.               C
!                                                                     C
!        YJ : Y COORDINATE OF THE POINT IN MODEL INDEX.               C
!                                                                     C
!-----------------------------------------------------------------------
 
      REAL                          :: XI, YJ, XLATP, XLONP
      REAL                          :: arg, arg1
      REAL                          :: deltalon, deltalon1
      REAL                          :: tl1r
      REAL                          :: rm
      REAL                          :: ctl1r
      REAL                          :: rebydx
      REAL                          :: cone, polei, polej, rsw
      integer                       :: i, j 

      rebydx = reby/dx*1e3

!     print *,' REBYDX : ',rebydx,dx

      ! Compute deltalon between known longitude and standard lon and ensure
      ! it is not in the cut zone
      deltalon = xlonp - stdlon
      IF (deltalon .GT. +180.) deltalon = deltalon - 360.
      IF (deltalon .LT. -180.) deltalon = deltalon + 360.
      deltalon1 = lon1 - stdlon
      IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
      IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.

!     print *,' DELTA : ',deltalon1, deltalon,lon1, stdlon

      ! Convert truelat1 to radian and compute COS for later use
      tl1r = truelat1 * rad_per_deg
      ctl1r = COS(tl1r)

      IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
         cone = ALOG10(COS(truelat1*rad_per_deg)) - &
                ALOG10(COS(truelat2*rad_per_deg))
         cone=cone/(ALOG10(TAN((45.0-ABS(truelat1)/2.0)*rad_per_deg))- &
                ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
      ELSE
         cone = SIN(ABS(truelat1)*rad_per_deg )
      ENDIF

!     print *,' tl1r : ',tl1r,truelat1, ctl1r,cone

      ! Radius to desired point
      rm = rebydx * ctl1r/cone * &
           (TAN((90.*hemi-xlatp)*rad_per_deg/2.) / &
            TAN((90.*hemi-truelat1)*rad_per_deg/2.))**cone

      arg = cone*(deltalon*rad_per_deg)
      arg1 = cone*(deltalon1*rad_per_deg)

      rsw = rebydx * ctl1r/cone * &
             (TAN((90.*hemi-lat1)*rad_per_deg/2.) / &
              TAN((90.*hemi-truelat1)*rad_per_deg/2.))**cone

!     print *,' RM : ',rm,arg,arg1,rsw

      polei = hemi*knowni - hemi * rsw * SIN(arg1)
      polej = hemi*knownj + rsw * COS(arg1)
      xi = polei + hemi * rm * SIN(arg)
      yj = polej - rm * COS(arg)

!     print *,' POLE : ',polei, polej, xi, yj, cos(arg1),cos(arg)

      ! Finally, if we are in the southern hemisphere, flip the i/j
      ! values to a coordinate system where (1,1) is the SW corner
      ! (what we assume) which is different than the original NCEP
      ! algorithms which used the NE corner as the origin in the
      ! southern hemisphere (left-hand vs. right-hand coordinate?)
      xi = hemi * xi
      yj = hemi * yj
      RETURN
      END SUBROUTINE  LAMBC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine check(status)
    integer, intent ( in) :: status
!   integer               :: nf_noerr

!   if(status /= nf_noerr) then
    if(status /= 0) then
!     print *, trim(nf_strerror(status))
      print *, 'NETCDF ERROR'
      stop "Stopped"
    end if
  end subroutine check

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      END MODULE MAP_GOCART_DMS_V01
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      PROGRAM GET_GOCART_BACKG_DMS

      USE MAP_GOCART_DMS_V01
      USE NETCDF

      CALL GET_GOCART_DMS_DATA

      END PROGRAM GET_GOCART_BACKG_DMS
