!
! $Header$
!
c
c
	subroutine read_pstoke(irec,
     .   zrec,zklono,zklevo,airefi,phisfi,
     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
     .   fm_therm,en_therm,
     .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)

C******************************************************************************
C  Frederic HOURDIN, Abderrahmane IDELKADI
C Lecture des parametres physique stockes online necessaires pour
C recalculer offline le transport de traceurs sur une grille 2x plus fine que 
C celle online 
C A FAIRE : une seule routine au lieu de 2 (lectflux, redecoupe)!
C******************************************************************************


       IMPLICIT NONE

#include "netcdf.inc"
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comgeom.h"
#include "temps.h"
#include "ener.h"
#include "logic.h"
#include "description.h"
#include "serre.h"
#include "indicesol.h"
#include "control.h"
#include "dimphy.h"
	
	  integer klono,klevo,imo,jmo
	  parameter (imo=iim/2,jmo=(jjm+1)/2)
	  parameter(klono=(jmo-1)*imo+2,klevo=llm)
	  REAL phisfi(klono)
          REAL phisfi2(imo,jmo+1),airefi2(imo,jmo+1)

          REAL mfu(klono,klevo), mfd(klono,klevo)
          REAL en_u(klono,klevo), de_u(klono,klevo)
          REAL en_d(klono,klevo), de_d(klono,klevo)
          REAL coefh(klono,klevo)
           REAL fm_therm(klono,klevo),en_therm(klono,klevo)

          REAL mfu2(imo,jmo+1,klevo), mfd2(imo,jmo+1,klevo)
          REAL en_u2(imo,jmo+1,klevo), de_u2(imo,jmo+1,klevo)
          REAL en_d2(imo,jmo+1,klevo), de_d2(imo,jmo+1,klevo)
          REAL coefh2(imo,jmo+1,klevo)
           REAL fm_therm2(imo,jmo+1,klevo)
           REAL en_therm2(imo,jmo+1,klevo)

          REAL pl(klevo)
          integer irec
          integer xid,yid,zid,tid
          real zrec,zklono,zklevo,zim,zjm
          integer ncrec,ncklono,ncklevo,ncim,ncjm

          real airefi(klono)
          character*20 namedim

c  !! attention !!
c attention il y a aussi le pb de def klono
c dim de phis??
	  
	 
          REAL frac_impa(klono,klevo), frac_nucl(klono,klevo)
          REAL frac_impa2(imo,jmo+1,klevo), 
     .     frac_nucl2(imo,jmo+1,klevo)
          REAL pyu1(klono), pyv1(klono)
          REAL pyu12(imo,jmo+1), pyv12(imo,jmo+1)
          REAL ftsol(klono,nbsrf)
          REAL psrf(klono,nbsrf) 
	  REAL ftsol1(klono),ftsol2(klono),ftsol3(klono),ftsol4(klono)
          REAL psrf1(klono),psrf2(klono),psrf3(klono),psrf4(klono)
          REAL ftsol12(imo,jmo+1),ftsol22(imo,jmo+1),
     .     ftsol32(imo,jmo+1),
     .     ftsol42(imo,jmo+1)
          REAL psrf12(imo,jmo+1),psrf22(imo,jmo+1),psrf32(imo,jmo+1),
     .     psrf42(imo,jmo+1)
		REAL t(klono,klevo)
		REAL t2(imo,jmo+1,klevo)	
	  integer ncidp
          save ncidp
		integer varidt
          integer varidmfu, varidmfd, varidps, varidenu, variddeu	
          integer varidend,varidded,varidch,varidfi,varidfn
           integer varidfmth,varidenth
          integer varidyu1,varidyv1,varidpl,varidai,varididvt
          integer varidfts1,varidfts2,varidfts3,varidfts4
          integer varidpsr1,varidpsr2,varidpsr3,varidpsr4
          save varidmfu, varidmfd, varidps, varidenu, variddeu
          save varidend,varidded,varidch,varidfi,varidfn
           save varidfmth,varidenth
          save varidyu1,varidyv1,varidpl,varidai,varididvt
          save varidfts1,varidfts2,varidfts3,varidfts4
          save varidpsr1,varidpsr2,varidpsr3,varidpsr4
		save varidt

          integer l, i
          integer start(4),count(4),status
          real rcode
          logical first
          save first
          data first/.true./



c ---------------------------------------------
c   Initialisation de la lecture des fichiers
c ---------------------------------------------

      if (irec .eq. 0) then

            ncidp=NCOPN('phystoke.nc',NCNOWRIT,rcode)

            varidps=NCVID(ncidp,'phis',rcode)
            print*,'ncidp,varidps',ncidp,varidps

            varidpl=NCVID(ncidp,'sig_s',rcode)
            print*,'ncidp,varidpl',ncidp,varidpl

            varidai=NCVID(ncidp,'aire',rcode)
            print*,'ncidp,varidai',ncidp,varidai

c A FAIRE: Es-il necessaire de stocke t?
	        varidt=NCVID(ncidp,'t',rcode)
                print*,'ncidp,varidt',ncidp,varidt

            varidmfu=NCVID(ncidp,'mfu',rcode)
            print*,'ncidp,varidmfu',ncidp,varidmfu

            varidmfd=NCVID(ncidp,'mfd',rcode)
            print*,'ncidp,varidmfd',ncidp,varidmfd

            varidenu=NCVID(ncidp,'en_u',rcode)
            print*,'ncidp,varidenu',ncidp,varidenu

            variddeu=NCVID(ncidp,'de_u',rcode)
            print*,'ncidp,variddeu',ncidp,variddeu

            varidend=NCVID(ncidp,'en_d',rcode)
            print*,'ncidp,varidend',ncidp,varidend
	
            varidded=NCVID(ncidp,'de_d',rcode)
            print*,'ncidp,varidded',ncidp,varidded
	
            varidch=NCVID(ncidp,'coefh',rcode)
            print*,'ncidp,varidch',ncidp,varidch
	
c abder (pour thermiques)
             varidfmth=NCVID(ncidp,'fm_th',rcode)
             print*,'ncidp,varidfmth',ncidp,varidfmth

             varidenth=NCVID(ncidp,'en_th',rcode)
             print*,'ncidp,varidenth',ncidp,varidenth

	    varidfi=NCVID(ncidp,'frac_impa',rcode)
            print*,'ncidp,varidfi',ncidp,varidfi
	
	    varidfn=NCVID(ncidp,'frac_nucl',rcode)
            print*,'ncidp,varidfn',ncidp,varidfn
	
            varidyu1=NCVID(ncidp,'pyu1',rcode)
            print*,'ncidp,varidyu1',ncidp,varidyu1
	
            varidyv1=NCVID(ncidp,'pyv1',rcode)
            print*,'ncidp,varidyv1',ncidp,varidyv1
	
            varidfts1=NCVID(ncidp,'ftsol1',rcode)
            print*,'ncidp,varidfts1',ncidp,varidfts1
	
            varidfts2=NCVID(ncidp,'ftsol2',rcode)
            print*,'ncidp,varidfts2',ncidp,varidfts2
         
            varidfts3=NCVID(ncidp,'ftsol3',rcode)
            print*,'ncidp,varidfts3',ncidp,varidfts3
  
            varidfts4=NCVID(ncidp,'ftsol4',rcode)
            print*,'ncidp,varidfts4',ncidp,varidfts4
	
            varidpsr1=NCVID(ncidp,'psrf1',rcode)
            print*,'ncidp,varidpsr1',ncidp,varidpsr1
	
            varidpsr2=NCVID(ncidp,'psrf2',rcode)
            print*,'ncidp,varidpsr2',ncidp,varidpsr2
	
	    varidpsr3=NCVID(ncidp,'psrf3',rcode)
            print*,'ncidp,varidpsr3',ncidp,varidpsr3

            varidpsr4=NCVID(ncidp,'psrf4',rcode)
            print*,'ncidp,varidpsr4',ncidp,varidpsr4
	
c ID pour les dimensions

            status = nf_inq_dimid(ncidp,'y',yid)
            status = nf_inq_dimid(ncidp,'x',xid)
            status = nf_inq_dimid(ncidp,'sig_s',zid)
            status = nf_inq_dimid(ncidp,'time_counter',tid)

c lecture des dimensions

            status = nf_inq_dim(ncidp,yid,namedim,ncjm)
            status = nf_inq_dim(ncidp,xid,namedim,ncim)
            status = nf_inq_dim(ncidp,zid,namedim,ncklevo)
            status = nf_inq_dim(ncidp,tid,namedim,ncrec)
	
            zrec=ncrec
            zklevo=ncklevo
            zim=ncim
            zjm=ncjm
	
	    zklono=zim*(zjm-2)+2
	
	    write(*,*) 'read_pstoke : zrec = ', zrec
            write(*,*) 'read_pstoke : zklevo = ', zklevo
            write(*,*) 'read_pstoke : zim = ', zim 
            write(*,*) 'read_pstoke : zjm = ', zjm
            write(*,*) 'read_pstoke : zklono = ', zklono

c niveaux de pression
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidpl,1,zklevo,pl)
#else
      status=NF_GET_VARA_REAL(ncidp,varidpl,1,zklevo,pl)
#endif

c lecture de aire et phis
	
      start(1)=1
      start(2)=1
      start(3)=1
      start(4)=0

      count(1)=zim
      count(2)=zjm
      count(3)=1
      count(4)=0

c phis
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidps,start,count,phisfi2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,phisfi2,phisfi)

c aire
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidai,start,count,airefi2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,airefi2,airefi)
      else

      print*,'ok1'

c ---------------------
c   lecture des champs
c ---------------------
	
	print*,'WARNING!!! Il n y a pas de test de coherence'
        print*,'sur le nombre de niveaux verticaux dans le fichier nc'

      start(1)=1
      start(2)=1
      start(3)=1
      start(4)=irec

      count(1)=zim
      count(2)=zjm
      count(3)=zklevo
      count(4)=1


C *** Lessivage******************************************************
c frac_impa 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfi,start,count,frac_impa2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfi,start,count,frac_impa2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_impa2,frac_impa)

c frac_nucl 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfn,start,count,frac_nucl2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfn,start,count,frac_nucl2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,frac_nucl2,frac_nucl)

C*** Temperature ******************************************************
c abder t
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidt,start,count,t2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,t2,t)

C*** Flux pour le calcul de la convection TIEDTK ***********************
c mfu
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidmfu,start,count,mfu2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfu2,mfu)

c mfd
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidmfd,start,count,mfd2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,mfd2,mfd)

c en_u 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidenu,start,count,en_u2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidenu,start,count,en_u2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_u2,en_u)

c de_u 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,variddeu,start,count,de_u2)
#else
      status=NF_GET_VARA_REAL(ncidp,variddeu,start,count,de_u2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_u2,de_u)

c en_d 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidend,start,count,en_d2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidend,start,count,en_d2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_d2,en_d)

c de_d 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidded,start,count,de_d2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidded,start,count,de_d2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,de_d2,de_d)

C **** Coeffecient du mellange turbulent**********************************
c coefh 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidch,start,count,coefh2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
#endif
       call gr_ecrit_fi(klevo,klono,imo,jmo+1,coefh2,coefh)

C*** Flux ascendant et entrant pour les Thermiques************************
cabder thermiques
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfmth,start,count,fm_therm2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfmth,start,count,fm_therm2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,fm_therm2,fm_therm)

#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidenth,start,count,en_therm2)
#else
      status=NF_GET_VARA_REAL(ncidp,varidenth,start,count,en_therm2)
#endif
      call gr_ecrit_fi(klevo,klono,imo,jmo+1,en_therm2,en_therm)

C*** Vitesses aux sol ******************************************************
      start(3)=irec
      start(4)=0
      count(3)=1
      count(4)=0
c pyu1
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidyu1,start,count,pyu12)
#else
      status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,pyu12,pyu1)

c pyv1
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidyv1,start,count,pyv12)
#else
      status=NF_GET_VARA_REAL(ncidp,varidyv1,start,count,pyv12)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,pyv12,pyv1)

C*** Temperature au sol ********************************************
c ftsol1
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfts1,start,count,ftsol12)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfts1,start,count,ftsol12)
#endif
       call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol12,ftsol1)

c ftsol2
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfts2,start,count,ftsol22)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfts2,start,count,ftsol22)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol22,ftsol2)

c ftsol3
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfts3,start,count,ftsol32)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfts3,start,count,ftsol32)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol32,ftsol3)

c ftsol4
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidfts4,start,count,ftsol42)
#else
      status=NF_GET_VARA_REAL(ncidp,varidfts4,start,count,ftsol42)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,ftsol42,ftsol4)

C*** Nature du sol **************************************************
c psrf1 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr1,start,count,psrf12)
#else
      status=NF_GET_VARA_REAL(ncidp,varidpsr1,start,count,psrf12)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf12,psrf1)

c psrf2 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr2,start,count,psrf22)
#else
      status=NF_GET_VARA_REAL(ncidp,varidpsr2,start,count,psrf22)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf22,psrf2)

c psrf3 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr3,start,count,psrf32)
#else
      status=NF_GET_VARA_REAL(ncidp,varidpsr3,start,count,psrf32)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf32,psrf3)

c psrf4 
#ifdef NC_DOUBLE
      status=NF_GET_VARA_DOUBLE(ncidp,varidpsr4,start,count,psrf42)
#else
      status=NF_GET_VARA_REAL(ncidp,varidpsr4,start,count,psrf42)
#endif
      call gr_ecrit_fi(1,klono,imo,jmo+1,psrf42,psrf4)
	
          do i = 1,klono
	
	psrf(i,1) = psrf1(i)
        psrf(i,2) = psrf2(i)
        psrf(i,3) = psrf3(i)
        psrf(i,4) = psrf4(i)
  
        ftsol(i,1) = ftsol1(i)
        ftsol(i,2) = ftsol2(i)
        ftsol(i,3) = ftsol3(i)
        ftsol(i,4) = ftsol4(i)
	
          enddo
	
	endif
	
	return
	
	end

