Index: /LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90
===================================================================
--- /LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 4356)
+++ /LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90	(revision 4357)
@@ -2,9 +2,12 @@
 ! $Id: $
 !
-! This subroutine creates the file grilles_gcm.nc containg longitudes and
-! latitudes in degrees for grid u and v. This subroutine is called from
-! ce0l. This subroutine corresponds to the first 
-! part in the program create_fausse_var.
+! 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)
 
@@ -20,32 +23,35 @@
   INCLUDE "netcdf.inc"
 
-
+!======================== 
   REAL,DIMENSION(iip1,jjp1),INTENT(IN)  :: masque ! masque terre/mer
   REAL,DIMENSION(iip1,jjp1),INTENT(IN)  :: phis   ! geopotentiel au sol
 
-  REAL temp(iim+1,jjm+1)
-  ! Attributs netcdf sortie
+  INTEGER status,i,j
+  
+  ! Attributs netcdf output
   INTEGER ncid_out,rcode_out
+  
   INTEGER out_lonuid,out_lonvid,out_latuid,out_latvid
-  INTEGER out_varid
+  INTEGER out_uid,out_vid,out_tempid
   INTEGER out_lonudim,out_lonvdim
   INTEGER out_latudim,out_latvdim,out_dim(2)
   INTEGER out_levdim
-
-  INTEGER start(4),COUNT(4)
-
-  INTEGER status,i,j
-  REAL rlatudeg(jjp1),rlatvdeg(jjm),rlevdeg(llm)
+  !
+  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 dlon1(iip1),dlon2(iip1),dlat1(jjp1),dlat2(jjp1)
-  REAL acoslat,dxkm,dykm,resol(iip1,jjp1)
-  REAL,DIMENSION(iip1,jjp1)  :: phis_loc
+  REAL uwnd(iip1,jjp1),vwnd(iip1,jjm),temp(iip1,jjp1) 
+  !
   INTEGER masque_int(iip1,jjp1)
-  INTEGER :: phis_id
-  INTEGER :: area_id
-  INTEGER :: mask_id
-  INTEGER :: presnivs_id
-  
+  REAL :: phis_loc(iip1,jjp1)
+  
+  !======================== 
+  ! CALCULATION of latu, latv, lonu, lonv in deg.
+  ! ---------------------------------------------------
   rad = 6400000
   omeg = 7.272205e-05
@@ -65,4 +71,5 @@
      rlatudeg(j)=rlatu(j)*180./pi
   ENDDO
+  
   DO j=1,jjm
      rlatvdeg(j)=rlatv(j)*180./pi
@@ -73,12 +80,31 @@
      rlonvdeg(i)=rlonv(i)*180./pi + 360.
   ENDDO
-
-
-  !  2 ----- OUVERTURE DE LA SORTIE NETCDF
-  ! ---------------------------------------------------
-  ! CREATION OUTPUT
-  ! ouverture fichier netcdf de sortie out
+ 
+  ! 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)
@@ -89,145 +115,118 @@
   status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)
   CALL handle_err(status)
-
-
-  !   Longitudes en u
+  !
+  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',9,'Longitude en u')
-
-  !   Longitudes en v
+  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', 9,'Longitude en v')
-
-  !   Latitude en u
+  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', 8,'Latitude en u')
-
-  !  Latitude en v
+  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', 8,'Latitude en v')
-
-  !   ecriture de la grille u
+  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_varid)
-  CALL handle_err(status)
-  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
-  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name', 16,'Grille aux point u')
-
-  !   ecriture de la grille v
+  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_varid)
-  CALL handle_err(status)
-  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
-  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name', 16,'Grille aux point v')
-
-  !   ecriture de la grille u
+  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_varid)
-  CALL handle_err(status)
-  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'units', 6,'Kelvin')
-  status=NF_PUT_ATT_TEXT(ncid_out,out_varid,'long_name',16,'Grille aux point u')
-
-  status=NF_ENDDEF(ncid_out)
-  write(*,*) "COUCOU 6"
-  CALL handle_err(status)
-  ! 5) ----- FERMETURE DES FICHIERS NETCDF------------------
-  ! --------------------------------------------------------
-  ! 3-b- Ecriture de la grille pour la sortie
-  ! rajoute l'ecriture de la grille
-
-  status=NF90_PUT_VAR(ncid_out,out_lonuid,rlonudeg,[1],[iim+1])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,out_lonvid,rlonvdeg,[1],[iim+1])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,out_latuid,rlatudeg,[1],[jjm+1])
-  CALL handle_err(status)
-  status=NF90_PUT_VAR(ncid_out,out_latvid,rlatvdeg,[1],[jjm])
-  CALL handle_err(status)
-
-  start(1)=1
-  start(2)=1
-  start(3)=1
-  start(4)=1
-
-  COUNT(1)=iim+1
-  COUNT(2)=jjm+1
-  COUNT(3)=1
-  COUNT(4)=1
-
-  DO j=1,jjm+1
-     DO i=1,iim+1
-        temp(i,j)=MOD(i,2)+MOD(j,2)
-     ENDDO
-  ENDDO
-
-  status=NF90_PUT_VAR(ncid_out,out_varid,temp,start, count)
-  CALL handle_err(status)
-
-  ! On re-ouvre le fichier pour rajouter 4 nouvelles variables necessaire pour INCA
-! lev - phis - aire - mask
-! rlevdeg(:) = presnivs
-  rlevdeg(:) = presnivs(:)
-  phis_loc(:,:) = phis(:,:)/g
-
-! niveaux de pression verticaux
-  status = NF_REDEF (ncid_out)
-  CALL handle_err(status)
-  status=NF_DEF_DIM(ncid_out,'lev',llm,out_levdim)
-  CALL handle_err(status)
-  status=NF90_DEF_VAR(ncid_out,'presnivs',NF90_FLOAT,out_levdim,&
-                    presnivs_id)
-  CALL handle_err(status)
-  
-! fields
+  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_ENDDEF(ncid_out)
-  CALL handle_err(status)
-
-  ! ecriture des variables
-  status=NF90_PUT_VAR(ncid_out,presnivs_id,rlevdeg,[1],[llm])
-  CALL handle_err(status)
-
-  start(1)=1
-  start(2)=1
-  start(3)=1
-  start(4)=0
+  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
-  COUNT(3)=1
-  COUNT(4)=0
-
+  
+  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)
+  CALL handle_err(status) 
   status = nf90_put_var(ncid_out, area_id, aire,start,count)
-  CALL handle_err(status)
-  masque_int(:,:) = nINT(masque(:,:))
+  CALL handle_err(status) 
   status = nf90_put_var(ncid_out, mask_id,masque_int,start,count)
   CALL handle_err(status)
   
-  ! fermeture du fichier netcdf
+  ! CLOSE netcdf file
   CALL ncclos(ncid_out,rcode_out)
+  write(*,*) "END grilles_gcm_netcdf_sub OK"
 
 END SUBROUTINE grilles_gcm_netcdf_sub
-
 
 
