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

      MODULE MAP_EROD_DATA_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
!
! ifort -O2 -FR -convert big_endian map_erod_v01.F -I /usr/local/netcdfINTEL/include -L/usr/local/netcdfINTEL/lib -lnetcdf
!
! compile with:
!  pgf90 -w -byteswapio -Mfree map_erod_v01.F
!
!or
!
!  mpif90 -axP -free -convert big_endian map_erod_v01.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 EROD DATA
!
!     PARAMETER( lon = 288 , ilon = 289 , lat = 181 , ilat = 182 , ndcls = 3 , ndsrc = 1 )
      integer       ::  lon = 288
      integer       :: ilon = 289
      integer       ::  lat = 181
      integer       :: ilat = 182
      integer       :: ndcls = 3
      integer       :: ndsrc = 1
      integer                        :: status,system


! - Output data dimensions for the WRF domain(user specified)
!
!     PARAMETER(IX2=39,JX2=39,KX=19,KP=KX+1)
      integer       ::  ix2 = 39
      integer       ::  jx2 = 39
      integer       ::  kx =  19
      integer       ::  kp =  20         ! kp = kx+1

!
      character (len=80)   ::  dir_in='/home/wrfchem/STEVEN/'
      character (len=80)   ::  fname_in_1='GAO_source_3cl.nc'
      character (len=80)   ::  fname_in_2='dust_asia.nc'

      character (len=80)   ::  fname_out='wrfchem_eroddata_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 = 60.E3
      real       :: lat1 =  38.00
      real       :: lon1 = -080.00
      real       :: stdlon = -80.0
      real       :: knowni = 19.0
      real       :: knownj = 19.0
      real       :: truelat1 = 38.0001
      real       :: truelat2 = 38.
      integer    :: hemi = +1

      CONTAINS

      SUBROUTINE  GET_EROD_DATA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      include 'netcdf.inc'

      real, dimension(lat)                     :: lat_g
      real, dimension(lon)                     :: lon_g
!     real, dimension(lat,lon)                 :: erod_mod_global, erod_total

!     real, dimension(ndsrc, ndcls, lat, lon)  :: erod_g
      real, dimension(lon, lat, ndcls, ndsrc )  :: erod_g
!     real, dimension(ndcls, lat, lon)  :: erod_g

      real, dimension(ix2,jx2,ndcls)         :: erod
!     real, dimension(ix2,jx2)         :: xlat,xlon
      real                             :: xlat,xlon
      real                             :: dist

      integer                        :: lat_id,lon_id

      integer     ::  i, j, k, jmx
      integer     ::  i2, j2
      real        ::  x2, y2

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

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

! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file.
      print *,' Open '
!     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) )

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

! Read the data.
      print *,' Lat  '
      call check( nf_get_var_real(ncid, lat_id, lat_g) )

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

! Read the data.
      print *,' Lon  '
      call check( nf_get_var_real(ncid, lon_id, lon_g) )

! Get the varid of the data variable, based on its name.
      print *,' Erod id  '
      call check( nf_inq_varid(ncid, "EROD", varid) )

! Read the data.
      print *,' Erod   '
      call check( nf_get_var_real(ncid, varid, erod_g) )

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Open the file. NF_NOWRITE tells netCDF we want read-only access to the file.
!      call check( nf_open(TRIM(dir_in)//TRIM(fname_in_2), NF_NOWRITE, ncid) )
!
!! Get the varid of the data variable, based on its name.
!      call check( nf_inq_varid(ncid, "EROD_MOD", varid) )
!
!! Read the data.
!      call check( nf_get_var(ncid, varid, erod_mod_global) )
!
!! Close the file, freeing all resources.
!      call check( nf_close(ncid) )

!  lat_g = real( lat_double)
!  lon_g = real( lon_double)

   print *,' LAT_G: ', lat_g(110:120)
   print * ,' '
   print *,' LON_G: ', lon_g(150:160)
   print * ,' '
!  print *,' EROD 1 ', erod_g(1,1,110:120,150:160)
   print *,' EROD 1 ', erod_g(150:160,110:120,1,1)
   print * ,' '

! apportion errosion map grid  onto wrf grid by simple grid dumping

!     real, dimension(lon, lat, ndcls, ndsrc )  :: erod_g

      do j=1,lat
      do i=1,lon

!  print * ,' CALL LAMBC', i, j ,' LAT, LON, ', lat_g(j), lon_g(i)

           CALL LAMBC(lat_g(j),lon_g(i),x2,y2)

!  print *,' LAT, LON, x, y:  ', lat_g(j), lon_g(i),x2, y2,i,j
!  print * ,' '


           if (X2>=1. .AND.  X2<=IX2 .AND.  Y2>=1. .AND.  Y2<=JX2 ) THEN
              i2=INT(x2)
              j2=INT(y2)

!  print *,' LATG, LONG:  ', lat_g(j), lon_g(i), x2, y2, erod_g(i,j,1,1)

               CALL MAPCF(float(i2),float(j2), xlat, xlon)
!  print *,' LAT , LON :  ', xlat, xlon, x2, y2
!  print * ,' '

              erod(i2, j2, 1) = erod_g(i, j, 1, 1)
              erod(i2, j2, 2) = erod_g(i, j, 2, 1)
              erod(i2, j2, 3) = erod_g(i, j, 3, 1)
            endif

       enddo
       enddo


! bogus some dust data for now

      do j=1,jx2
      do i=1,ix2
         dist = (float(i) - float(ix2/2)) * (float(i) - float(ix2/2)) + &
                (float(j) - float(jx2/2)) * (float(j) - float(jx2/2)) 

         erod(i,j,1) = 1.E+9 * exp(-1*dist/3.)
         erod(i,j,2) = 2.E+9 * exp(-1*dist/3.)
         erod(i,j,3) = 3.E+9 * exp(-1*dist/3.)

         if(erod(i,j,1) < 1.e5) erod(i,j,1) = 0.
         if(erod(i,j,2) < 1.e5) erod(i,j,2) = 0.
         if(erod(i,j,3) < 1.e5) erod(i,j,3) = 0.
      enddo
      enddo

   print *,' EROD  ', erod(1:40,20:20,1)
   print * ,' '

! Open output files
      OPEN(19,FILE='wrf_dust_erod.txt',FORM='FORMATTED')
      write(19,'(12E15.5)') erod(:,:,1)
      write(19,'(12E15.5)') erod(:,:,2)
      write(19,'(12E15.5)') erod(:,:,3)

      CLOSE(19)
      END SUBROUTINE  GET_EROD_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
      INTEGER :: jnew, inew
      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 * i
      jnew = hemi * j

      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.

      write(*,*) 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

      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

    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_EROD_DATA_V01
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      PROGRAM GET_ERROSION_DATA

      USE MAP_EROD_DATA_V01
      USE NETCDF


      CALL GET_EROD_DATA

      END PROGRAM GET_ERROSION_DATA
