!
! $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
  
  IMPLICIT NONE

  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"
  INCLUDE "comgeom.h"
  INCLUDE "netcdf.inc"

!======================== 
  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=NF_CREATE('grilles_gcm.nc',IOR(NF_CLOBBER,NF_64BIT_OFFSET),ncid_out)
  CALL handle_err(status)
  
  ! DEFINE output dimensions
  !-------------------------
  status=NF_DEF_DIM(ncid_out,'lonu',iim+1,out_lonudim)
  CALL handle_err(status)
  status=NF_DEF_DIM(ncid_out,'lonv',iim+1,out_lonvdim)
  CALL handle_err(status)
  status=NF_DEF_DIM(ncid_out,'latu',jjm+1,out_latudim)
  CALL handle_err(status)
  status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)
  CALL handle_err(status)
  !
  status=NF_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=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'units', 12,'degrees_east')
  status=NF_PUT_ATT_TEXT(ncid_out,out_lonuid,'long_name',19,'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=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'units', 12,'degrees_east')
  status=NF_PUT_ATT_TEXT(ncid_out,out_lonvid,'long_name', 19,'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=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'units', 13,'degrees_north')
  status=NF_PUT_ATT_TEXT(ncid_out,out_latuid,'long_name', 18,'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=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'units', 13,'degrees_north')
  status=NF_PUT_ATT_TEXT(ncid_out,out_latvid,'long_name', 18,'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=NF_PUT_ATT_TEXT(ncid_out,out_uid,'units', 3,'m/s')
  status=NF_PUT_ATT_TEXT(ncid_out,out_uid,'long_name', 21,'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=NF_PUT_ATT_TEXT(ncid_out,out_vid,'units', 3,'m/s')
  status=NF_PUT_ATT_TEXT(ncid_out,out_vid,'long_name', 21,'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=NF_PUT_ATT_TEXT(ncid_out,out_tempid,'units', 6,'Kelvin')
  status=NF_PUT_ATT_TEXT(ncid_out,out_tempid,'long_name',21,'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=NF_PUT_ATT_TEXT(ncid_out,presnivs_id,'units', 2,'Pa')
  status=NF_PUT_ATT_TEXT(ncid_out,presnivs_id,'long_name',15,'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=NF_PUT_ATT_TEXT(ncid_out,phis_id,'units', 1,'m')
  status=NF_PUT_ATT_TEXT(ncid_out,phis_id,'long_name',27,'surface geopotential height')
  ! gridcell area
  status = nf90_def_var(ncid_out,'aire',NF90_FLOAT,out_dim,area_id)
  CALL handle_err(status)
  status=NF_PUT_ATT_TEXT(ncid_out,area_id,'units', 2,'m2')
  status=NF_PUT_ATT_TEXT(ncid_out,area_id,'long_name',13,'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=NF_PUT_ATT_TEXT(ncid_out,mask_id,'long_name',27,'land-sea mask (nINT approx)')
  
  ! END the 'define' mode in netCDF file
  status=NF_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)
  INCLUDE "netcdf.inc"

  INTEGER status
  IF (status.NE.nf_noerr) THEN
     PRINT *,NF_STRERROR(status)
     CALL abort_gcm('grilles_gcm_netcdf','netcdf error',1)
  ENDIF
END SUBROUTINE handle_err