!
! $Id: grilles_gcm_netcdf.F 1907 2013-11-26 13:10:46Z emillour $
!
c
c

      PROGRAM create_fausse_var
C
      IMPLICIT NONE
C
C
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comgeom.h"
#include "comvert.h"

      real temp(iim+1,jjm+1)
#include "netcdf.inc"

c Attributs netcdf sortie
        character*64 fich_out
        integer*4 ncid_out,rcode_out
        integer*4 out_lonuid,out_lonvid,out_latuid,out_latvid
        integer*4 out_varid
        integer*4 out_lonudim,out_lonvdim
        integer*4 out_latudim,out_latvdim,out_dim(3)

      INTEGER         longcles
      PARAMETER     ( longcles = 20 )
      REAL  clesphy0( longcles )

      integer start(4),count(4)

	integer status,i,j
        real rlatudeg(jjp1),rlatvdeg(jjm)
        real rlonudeg(iip1),rlonvdeg(iip1)

      real dlon1(iip1),dlon2(iip1),dlat1(jjp1),dlat2(jjp1)
      real acoslat,dxkm,dykm,resol(iip1,jjp1)

#include "serre.h"
#include "fxyprim.h"

      print*,'OK0'

      rad = 6400000
      omeg = 7.272205e-05
      g = 9.8
      kappa = 0.285716
      daysec = 86400
      cpp = 1004.70885

      preff = 101325.
      pa= 50000.

c     open(99,file='run.def',status='old',form='formatted')
c     CALL defrun_new( 99, .TRUE.,clesphy0 )
c     close(99)

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


      print*,'OK1'
      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


      print*,'OK2'
c  2 ----- OUVERTURE DE LA SORTIE NETCDF
c ---------------------------------------------------
c CREATION OUTPUT
c ouverture fichier netcdf de sortie out
        fich_out='grilles_gcm.nc'

        status=NF_CREATE(fich_out,NF_NOCLOBBER,ncid_out)
        status=NF_DEF_DIM(ncid_out,'lonu',iim+1,out_lonudim)
        status=NF_DEF_DIM(ncid_out,'lonv',iim+1,out_lonvdim)
        status=NF_DEF_DIM(ncid_out,'latu',jjm+1,out_latudim)
        status=NF_DEF_DIM(ncid_out,'latv',jjm,out_latvdim)


      print*,'OK3'
c   Longitudes en u
        print *,'OUTID: ',ncid_out
        status=NF_DEF_VAR(ncid_out,'lonu',NF_FLOAT,1,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')

c   Longitudes en v
        print *,'OUTID: ',ncid_out
        status=NF_DEF_VAR(ncid_out,'lonv',NF_FLOAT,1,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')

c   Latitude en u
        status=NF_DEF_VAR(ncid_out,'latu',NF_FLOAT,1,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')

c  Latitude en v
        status=NF_DEF_VAR(ncid_out,'latv',NF_FLOAT,1,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')

c   ecriture de la grille u
        out_dim(1)=out_lonudim
        out_dim(2)=out_latudim
        status=NF_DEF_VAR(ncid_out,'grille_u',NF_FLOAT,2,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')

c   ecriture de la grille v
        out_dim(1)=out_lonvdim
        out_dim(2)=out_latvdim
        status=NF_DEF_VAR(ncid_out,'grille_v',NF_FLOAT,2,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')

c   ecriture de la grille u
        out_dim(1)=out_lonvdim
        out_dim(2)=out_latudim
        status=NF_DEF_VAR(ncid_out,'grille_s',NF_FLOAT,2,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')


      print*,'OK4'
        status=NF_ENDDEF(ncid_out)
c 5) ----- FERMETURE DES FICHIERS NETCDF------------------
c --------------------------------------------------------
c 3-b- Ecriture de la grille pour la sortie
c rajoute l'ecriture de la grille

#ifdef NC_DOUBLE
      status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonuid,1,iim+1,rlonudeg)
      status=NF_PUT_VARA_DOUBLE(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
      status=NF_PUT_VARA_DOUBLE(ncid_out,out_latuid,1,jjm+1,rlatudeg)
      status=NF_PUT_VARA_DOUBLE(ncid_out,out_latvid,1,jjm,rlatvdeg)
#else
      status=NF_PUT_VARA_REAL(ncid_out,out_lonuid,1,iim+1,rlonudeg)
      status=NF_PUT_VARA_REAL(ncid_out,out_lonvid,1,iim+1,rlonvdeg)
      status=NF_PUT_VARA_REAL(ncid_out,out_latuid,1,jjm+1,rlatudeg)
      status=NF_PUT_VARA_REAL(ncid_out,out_latvid,1,jjm,rlatvdeg)
#endif

        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

#ifdef NC_DOUBLE
        status=NF_PUT_VARA_DOUBLE(ncid_out,out_varid,start,
     s  count,temp)
#else
        status=NF_PUT_VARA_REAL(ncid_out,out_varid,start,
     s  count,temp)
#endif


c fermeture du fichier netcdf
        call ncclos(ncid_out,rcode_out)
        write(*,*) 'Fermeture: ',fich_out


      print*,'OK5'
c   Ecriture grads
      open (20,file='grille.dat',form='unformatted',access='direct'
     s      ,recl=4*ip1jmp1)
      write(20,rec=1) ((REAL(mod(i,2)+mod(j,2)),i=1,iip1),j=1,jjp1)
      write(20,rec=2) ((REAL(mod(i,2)*mod(j,2)),i=1,iip1),j=1,jjp1)
      do j=2,jjm
         dlat1(j)=180.*(rlatv(j)-rlatv(j-1))/pi
c        dlat2(j)=180.*fyprim(REAL(j))/pi
      enddo
      do i=2,iip1
         dlon1(i)=180.*(rlonu(i)-rlonu(i-1))/pi
c        dlon2(i)=180.*fxprim(REAL(i))/pi
      enddo
      do j=2,jjm
         dykm=(rlatv(j)-rlatv(j-1))*6400.
         acoslat=6400.*cos(rlatu(j))
         do i=2,iip1
            dxkm=acoslat*(rlonu(i)-rlonu(i-1))
            resol(i,j)=sqrt(dykm*dykm+dxkm*dxkm)
         enddo
         resol(1,j)=resol(iip1,j)
      enddo
      write(20,rec=3) resol
      dlon1(1)=dlon1(iip1)
      dlon2(1)=dlon2(iip1)
      write(20,rec=4) ((dlon1(i),i=1,iip1),j=1,jjp1)
      write(20,rec=5) ((dlon1(i)*pi/180.*0.001*
     s   cos(rlatu(j))*rad,i=1,iip1),j=1,jjp1)
      write(20,rec=6) ((dlon2(i),i=1,iip1),j=1,jjp1)
      write(20,rec=7) ((dlat1(j),i=1,iip1),j=1,jjp1)
      write(20,rec=8) ((dlat1(j)*pi/180.*rad*0.001,i=1,iip1),j=1,jjp1)
      write(20,rec=9) ((dlat2(j),i=1,iip1),j=1,jjp1)

      print*,'I, LON, DX (km)'
      do i=1,iip1
         print*,i,rlonu(i)*180./pi,dlon1(i)*pi/180.*0.001*
     s   cos(clat*pi/180.)*rad
      enddo
      print*,'J, LAT, DY (km)'
      do j=1,jjp1
         print*,j,rlatu(j)*180./pi,dlat1(j)*pi/180.*0.001*rad
      enddo

      open (21,file='grille.ctl',form='formatted')

c   WARNING! on reecrase le fichier .ctl a chaque ecriture
      write(21,'(a5,1x,a40)')
     &       'DSET ','^grille.dat'

      write(21,'(a12)') 'UNDEF 1.0E30'
      write(21,'(a5,1x,a40)') 'TITLE ','grille'
      call formcoord(21,iip1,rlonv,180./pi,.false.,'XDEF')
      call formcoord(21,jjp1,rlatu,180./pi,.true.,'YDEF')
      call formcoord(21,1,0.,1.,.false.,'ZDEF')
      write(21,'(a4,i10,a30)')
     &       'TDEF ',1,' LINEAR 23OCT1994 3hr '
      write(21,'(a4,2x,i5)') 'VARS',9
      write(21,'(a18)') 'grille 0 99 grille'
      write(21,'(a18)') 'gril   0 99 gril  '
      write(21,'(a29)') 'resol   0 99 resolution (km)  '
      write(21,'(a18)') 'dlon1  0 99 dlon1 '
      write(21,'(a20)') 'dx     0 99 dx (km) '
      write(21,'(a18)') 'dlon2  0 99 dlon2 '
      write(21,'(a18)') 'dlat1  0 99 dlat1 '
      write(21,'(a20)') 'dy     0 99 dy (km) '
      write(21,'(a18)') 'dlat2  0 99 dlat2 '
      write(21,'(a7)') 'ENDVARS'





      print*,'OK6'
	end



        subroutine handle_err(status)
#include "netcdf.inc"


        integer status
        print *,'handle code err: ',NF_NOERR
        IF (status.NE.nf_noerr) THEN
                print *,NF_STRERROR(status)
                stop 'stopped'
        ENDIF
        END