c
c $Header$
c
      SUBROUTINE lectfluxnc(irec,massem,pbarun,pbarvn,wn,tetan,phin,
     s     nrec,avant,airefi,phisfi,
     s     t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
     s     yu1,yv1,ftsol,pctsrf,
     s     frac_impa,frac_nucl,phis)

      IMPLICIT NONE

#include "dimensions.h"
#include "paramet.h"
#include "logic.h"
#include "comvert.h"
#include "comconst.h"
#include "comgeom2.h"

#include "tracstoke.h"

      integer irec,nrec,i,j

      integer klon,ig,l
      parameter (klon=iim*(jjm-1)+2)
      INTEGER nbsrf
      PARAMETER (nbsrf=4) ! nombre de sous-fractions pour une maille

c Convection
      real*4 zmfd(klon,llm),zde_d(klon,llm),zen_d(klon,llm)
      real*4 zmfu(klon,llm),zde_u(klon,llm),zen_u(klon,llm)
      real*4 mfd(klon,llm),de_d(klon,llm),en_d(klon,llm)
      real*4 mfu(klon,llm),de_u(klon,llm),en_u(klon,llm)
	real*4 t(klon,llm)
      real zdtvr
	
	real*4 airedy(iip1,jjp1)
      real*4 rlonu_dy(iip1,jjp1),rlonv_dy(iip1,jjm),
     . rlatu_dy(iip1,jjp1),rlatv_dy(iip1,jjm)

      real*4 coefkz(klon,llm)
      real*4 frac_impa(klon,llm),frac_nucl(klon,llm)
      real*4 yu1(klon), yv1(klon)
      real*4 ftsol(klon,nbsrf),pctsrf(klon,nbsrf)

      logical avant

      real massefi(klon,llm)

c   Flux masse
      real*4 massem(iip1,jjp1,llm),tetan(iip1,jjp1,llm)
      real*4 pbarun(iip1,jjp1,llm),pbarvn(iip1,jjm,llm)
      real*4 pbarvst(iip1,jjp1,llm)
      real*4 wn(iip1,jjp1,llm),phin(iip1,jjp1,llm)
      real*4 phis(iip1,jjp1)

      real*4 airefi(klon),phisfi(klon)

      real zcontrole(klon),zmass,tmpdyn(iip1,jjp1),zflux

      real ziadvtrac,ziadvtrac2,zrec2
      integer zim,zjm,zlm,zklon,zklev,zrec

      real zpi

      zpi=2.*asin(1.)


c==================================================================
c   Si le numero du record est 0 alors: INITIALISATION
c==================================================================
c
      print*,'ENTREE DANS LECTFLUXNC'
        print*,'IREC=',IREC

      if(irec.eq.0) then

        print*,'IREC==',0

C test         call inigeom

c==================================================================
c   ouverture des fichiers
c==================================================================
	
	CALL read_dstoke(0,zdtvr,ziadvtrac,ziadvtrac2)
	
	CALL read_fstoke0(0,
     .  zrec,zim,zjm,zlm,
     .   rlonu_dy,rlonv_dy,rlatu_dy,rlatv_dy,airedy,phis,
     .   massem,pbarun,pbarvn,wn,tetan,phin)
c 
	if(physic)then	
	CALL read_pstoke0(0,
     .   zrec,zklon,zklev,airefi,phisfi,
     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefkz,
     .   frac_impa,frac_nucl,yu1,yv1,ftsol,pctsrf)
	endif
	
	dtvr=zdtvr	
	nrec=zrec
        print*,'nrec=',nrec
	
	print*, 'Lecture de defstoke.nc'
	print*, 'Ds lectfluxnc dtvr = ', dtvr
	print*, 'Ds lectfluxnc istdyn= ',ziadvtrac
        print*, 'Ds lectfluxnc istphy= ',ziadvtrac2

      istdyn=ziadvtrac
      istphy=ziadvtrac2
	
c	pause

c==================================================================
c   Fin des initialisations
      else ! irec=0
c==================================================================


c-----------------------------------------------------------------------
c   Lecture des fichiers fluxmass et  physique:
c   -----------------------------------------------------

	 CALL read_fstoke0(irec,
     .   zrec,zim,zjm,zlm,
     .   rlonu_dy,rlonv_dy,rlatu_dy,rlatv_dy,airedy,phis,
     .   massem,pbarun,pbarvn,wn,tetan,phin)
	
	write(*,*) 'lectfluxnc aps read_fstoke irec',irec
	if(physic)then
        CALL read_pstoke0(irec,
     .   zrec,zklon,zklev,airefi,phisfi,
     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefkz,
     .   frac_impa,frac_nucl,yu1,yv1,ftsol,pctsrf)
	
	write(*,*) 'lectfluxnc aps read_pstoke irec',irec
	
	if (avant) then
c Simu directe
       do l=1,llm
          do ig=1,klon
             zmfu(ig,l)=mfu(ig,l)
             zmfd(ig,l)=mfd(ig,l)
             zde_u(ig,l)=de_u(ig,l)
             zen_u(ig,l)=en_u(ig,l)
             zde_d(ig,l)=de_d(ig,l)
             zen_d(ig,l)=en_d(ig,l)
          enddo
       enddo
      else
c   Simu retro
       do l=1,llm
          do ig=1,klon
             zmfd(ig,l)=-mfu(ig,l)
             zmfu(ig,l)=-mfd(ig,l)
             zen_d(ig,l)=de_u(ig,l)
             zde_d(ig,l)=en_u(ig,l)
             zen_u(ig,l)=de_d(ig,l)
             zde_u(ig,l)=en_d(ig,l)
          enddo
       enddo
      endif

c-----------------------------------------------------------------------
c   PETIT CONTROLE SUR LES FLUX CONVECTIFS...
c-----------------------------------------------------------------------

      print*,'Ap redec irec'
	write(*,*) 'dtphys = ', dtphys

         call gr_dyn_fi(llm,iip1,jjp1,klon,massem,massefi)

         do ig=1,klon
            zcontrole(ig)=1.
         enddo
c   zmass=(max(massem(ig,l),massem(ig,l-1))/airefi(ig)
         do l=2,llm
            do ig=1,klon
               zmass=max(massefi(ig,l),massefi(ig,l-1))/airefi(ig)
               zflux=max(abs(zmfu(ig,l)),abs(zmfd(ig,l)))*dtphys
               if(zflux.gt.0.9*zmass) then
                 zcontrole(ig)=min(zcontrole(ig),0.9*zmass/zflux)
               endif
            enddo
         enddo

         do ig=1,klon
            if(zcontrole(ig).lt.0.99999) then
               print*,'ATTENTION !!! on reduit les flux de masse '
               print*,'convectifs au point ig=',ig
            endif
         enddo
	
	write(*,*) 'lectfluxnc avt gr_fi_dyn'

         call gr_fi_dyn(1,klon,iip1,jjp1,zcontrole,tmpdyn)
	
	write(*,*) 'lectfluxnc aps gr_fi_dyn'
        
         do l=1,llm
            do ig=1,klon
               zmfu(ig,l)=zmfu(ig,l)*zcontrole(ig)
               zmfd(ig,l)=zmfd(ig,l)*zcontrole(ig)
               zen_u(ig,l)=zen_u(ig,l)*zcontrole(ig)
               zde_u(ig,l)=zde_u(ig,l)*zcontrole(ig)
               zen_d(ig,l)=zen_d(ig,l)*zcontrole(ig)
               zde_d(ig,l)=zde_d(ig,l)*zcontrole(ig)
            enddo
         enddo
	 endif
	
	write(*,*) 'lectfluxnc fin de boucle if'

      endif ! irec=0

      RETURN
      END