! $Id: $

! This SUBROUTINE creates the grilles_gcm.nc file, containing:
! -> longitudes and latitudes in degrees for dynamical grids u, v and scalaire,
! and the following variables added for INCA (informative anyway)
! -> vertical levels "presnivs"
! -> mask (land/sea), area (grid), phis=surface geopotential height = phis/g

! The SUBROUTINE is called in dynphy_lonlat/phylmd/ce0l.F90.

SUBROUTINE grilles_gcm_netcdf_sub(masque, phis)

  USE comconst_mod, ONLY: cpp, kappa, g, omeg, daysec, rad, pi
  USE comvert_mod, ONLY: presnivs, preff, pa
  USE netcdf, ONLY: nf90_def_var, nf90_int, nf90_float, nf90_put_var, nf90_enddef, &
          nf90_put_att, nf90_def_dim, nf90_64bit_offset, nf90_clobber, nf90_create
  USE lmdz_comgeom

  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
  USE lmdz_paramet
  USE lmdz_conf_gcm, ONLY: conf_gcm

  IMPLICIT NONE


  !========================
  REAL, DIMENSION(iip1, jjp1), INTENT(IN) :: masque ! masque terre/mer
  REAL, DIMENSION(iip1, jjp1), INTENT(IN) :: phis   ! geopotentiel au sol

  INTEGER status, i, j

  ! Attributs netcdf output
  INTEGER ncid_out, rcode_out

  INTEGER out_lonuid, out_lonvid, out_latuid, out_latvid
  INTEGER out_uid, out_vid, out_tempid
  INTEGER out_lonudim, out_lonvdim
  INTEGER out_latudim, out_latvdim, out_dim(2)
  INTEGER out_levdim

  INTEGER :: presnivs_id
  INTEGER :: mask_id, area_id, phis_id

  INTEGER start(2), COUNT(2)

  ! Variables
  REAL rlatudeg(jjp1), rlatvdeg(jjm), rlev(llm)
  REAL rlonudeg(iip1), rlonvdeg(iip1)
  REAL uwnd(iip1, jjp1), vwnd(iip1, jjm), temp(iip1, jjp1)

  INTEGER masque_int(iip1, jjp1)
  REAL :: phis_loc(iip1, jjp1)

  !======================== 
  ! CALCULATION of latu, latv, lonu, lonv in deg.
  ! ---------------------------------------------------
  rad = 6400000
  omeg = 7.272205e-05
  g = 9.8
  kappa = 0.285716
  daysec = 86400
  cpp = 1004.70885

  preff = 101325.
  pa = 50000.

  CALL conf_gcm(99, .TRUE.)
  CALL iniconst
  CALL inigeom

  DO j = 1, jjp1
    rlatudeg(j) = rlatu(j) * 180. / pi
  ENDDO

  DO j = 1, jjm
    rlatvdeg(j) = rlatv(j) * 180. / pi
  ENDDO

  DO i = 1, iip1
    rlonudeg(i) = rlonu(i) * 180. / pi + 360.
    rlonvdeg(i) = rlonv(i) * 180. / pi + 360.
  ENDDO

  ! CALCULATION of "false" variables on u, v, s grids
  ! --------------------------------------------------- 
  DO i = 1, iip1
    DO j = 1, jjp1
      uwnd(i, j) = MOD(i, 2) + MOD(j, 2)
      temp(i, j) = MOD(i, 2) + MOD(j, 2)
    ENDDO
    DO j = 1, jjm
      vwnd(i, j) = MOD(i, 2) + MOD(j, 2)
    END DO
  ENDDO

  ! CALCULATION of local vars for presnivs, mask, sfc. geopot. height
  ! --------------------------------------------------- 
  rlev(:) = presnivs(:)
  phis_loc(:, :) = phis(:, :) / g
  masque_int(:, :) = nINT(masque(:, :))


  ! OPEN output netcdf file
  !-------------------------
  status = nf90_create('grilles_gcm.nc', IOR(nf90_clobber, nf90_64bit_offset), ncid_out)
  CALL handle_err(status)

  ! DEFINE output dimensions
  !-------------------------
  status = nf90_def_dim(ncid_out, 'lonu', iim + 1, out_lonudim)
  CALL handle_err(status)
  status = nf90_def_dim(ncid_out, 'lonv', iim + 1, out_lonvdim)
  CALL handle_err(status)
  status = nf90_def_dim(ncid_out, 'latu', jjm + 1, out_latudim)
  CALL handle_err(status)
  status = nf90_def_dim(ncid_out, 'latv', jjm, out_latvdim)
  CALL handle_err(status)

  status = nf90_def_dim(ncid_out, 'lev', llm, out_levdim)
  CALL handle_err(status)

  ! DEFINE output variables
  !-------------------------
  !   Longitudes on "u" dynamical grid
  status = nf90_def_var(ncid_out, 'lonu', nf90_float, out_lonudim, out_lonuid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_lonuid, 'units', 'degrees_east')
  status = nf90_put_att(ncid_out, out_lonuid, 'long_name', 'Longitude on u grid')
  !   Longitudes on "v" dynamical grid
  status = nf90_def_var(ncid_out, 'lonv', nf90_float, out_lonvdim, out_lonvid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_lonvid, 'units', 'degrees_east')
  status = nf90_put_att(ncid_out, out_lonvid, 'long_name', 'Longitude on v grid')
  !   Latitudes on "u" dynamical grid
  status = nf90_def_var(ncid_out, 'latu', nf90_float, out_latudim, out_latuid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_latuid, 'units', 'degrees_north')
  status = nf90_put_att(ncid_out, out_latuid, 'long_name', 'Latitude on u grid')
  !  Latitudes on "v" dynamical grid
  status = nf90_def_var(ncid_out, 'latv', nf90_float, out_latvdim, out_latvid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_latvid, 'units', 'degrees_north')
  status = nf90_put_att(ncid_out, out_latvid, 'long_name', 'Latitude on v grid')
  !  "u" lat/lon dynamical grid
  out_dim(1) = out_lonudim
  out_dim(2) = out_latudim
  status = nf90_def_var(ncid_out, 'grille_u', nf90_float, out_dim, out_uid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_uid, 'units', 'm/s')
  status = nf90_put_att(ncid_out, out_uid, 'long_name', 'u-wind dynamical grid')
  !  "v" lat/lon dynamical grid
  out_dim(1) = out_lonvdim
  out_dim(2) = out_latvdim
  status = nf90_def_var(ncid_out, 'grille_v', nf90_float, out_dim, out_vid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_vid, 'units', 'm/s')
  status = nf90_put_att(ncid_out, out_vid, 'long_name', 'v-wind dynamical grid')
  !  "s" (scalar) lat/lon dynamical grid
  out_dim(1) = out_lonvdim
  out_dim(2) = out_latudim
  status = nf90_def_var(ncid_out, 'grille_s', nf90_float, out_dim, out_tempid)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, out_tempid, 'units', 'Kelvin')
  status = nf90_put_att(ncid_out, out_tempid, 'long_name', 'scalar dynamical grid')

  ! for INCA :
  ! vertical levels "presnivs"
  status = nf90_def_var(ncid_out, 'presnivs', nf90_float, out_levdim, presnivs_id)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, presnivs_id, 'units', 'Pa')
  status = nf90_put_att(ncid_out, presnivs_id, 'long_name', 'Vertical levels')
  ! surface geopotential height: named "phis" as the sfc geopotential, is actually phis/g
  out_dim(1) = out_lonvdim
  out_dim(2) = out_latudim
  status = nf90_def_var(ncid_out, 'phis', nf90_float, out_dim, phis_id)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, phis_id, 'units', 'm')
  status = nf90_put_att(ncid_out, phis_id, 'long_name', 'surface geopotential height')
  ! gridcell area
  status = nf90_def_var(ncid_out, 'aire', nf90_float, out_dim, area_id)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, area_id, 'units', 'm2')
  status = nf90_put_att(ncid_out, area_id, 'long_name', 'gridcell area')
  ! land-sea mask (nearest integer approx)
  status = nf90_def_var(ncid_out, 'mask', nf90_int, out_dim, mask_id)
  CALL handle_err(status)
  status = nf90_put_att(ncid_out, mask_id, 'long_name', 'land-sea mask (nINT approx)')

  ! END the 'define' mode in netCDF file
  status = nf90_enddef(ncid_out)
  CALL handle_err(status)

  ! WRITE the variables
  !-------------------------
  ! 1D : lonu, lonv,latu,latv ; INCA : presnivs
  status = nf90_put_var(ncid_out, out_lonuid, rlonudeg, [1], [iip1])
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, out_lonvid, rlonvdeg, [1], [iip1])
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, out_latuid, rlatudeg, [1], [jjp1])
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, out_latvid, rlatvdeg, [1], [jjm])
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, presnivs_id, rlev, [1], [llm])
  CALL handle_err(status)

  ! 2D : grille_u,grille_v,grille_s ; INCA: phis,aire,mask
  start(:) = 1
  COUNT(1) = iip1

  COUNT(2) = jjp1  ! for "u" and "s" grids
  status = nf90_put_var(ncid_out, out_uid, uwnd, start, count)
  CALL handle_err(status)
  COUNT(2) = jjm  ! for "v" grid
  status = nf90_put_var(ncid_out, out_vid, vwnd, start, count)
  CALL handle_err(status)
  COUNT(2) = jjp1  ! as "s" grid, for all the following vars
  status = nf90_put_var(ncid_out, out_tempid, temp, start, count)
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, phis_id, phis_loc, start, count)
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, area_id, aire, start, count)
  CALL handle_err(status)
  status = nf90_put_var(ncid_out, mask_id, masque_int, start, count)
  CALL handle_err(status)

  ! CLOSE netcdf file
  CALL ncclos(ncid_out, rcode_out)
  WRITE(*, *) "END grilles_gcm_netcdf_sub OK"

END SUBROUTINE grilles_gcm_netcdf_sub


SUBROUTINE handle_err(status)
  USE netcdf, ONLY: nf90_strerror, nf90_noerr

  INTEGER status
  IF (status/=nf90_noerr) THEN
    PRINT *, nf90_strerror(status)
    CALL abort_gcm('grilles_gcm_netcdf', 'netcdf error', 1)
  ENDIF
END SUBROUTINE handle_err