Index: /LMDZ.3.3/trunk/libf/dyn3d/lectfluxnc.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/lectfluxnc.F	(revision 185)
+++ /LMDZ.3.3/trunk/libf/dyn3d/lectfluxnc.F	(revision 185)
@@ -0,0 +1,215 @@
+      SUBROUTINE lectfluxnc0(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
+
+
Index: /LMDZ.3.3/trunk/libf/dyn3d/read_pstoke.F
===================================================================
--- /LMDZ.3.3/trunk/libf/dyn3d/read_pstoke.F	(revision 185)
+++ /LMDZ.3.3/trunk/libf/dyn3d/read_pstoke.F	(revision 185)
@@ -0,0 +1,355 @@
+	subroutine read_pstoke0(irec,
+     .   zrec,zkon,zkev,airefi,phisfi,
+     .   t,mfu,mfd,en_u,de_u,en_d,de_d,coefh,
+     .   frac_impa,frac_nucl,pyu1,pyv1,ftsol,psrf)
+
+
+       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*4 kon,kev,zkon,zkev
+	  parameter(kon=iim*(jjm-1)+2,kev=llm)
+	  REAL*4 phisfi(kon)
+          REAL*4 phisfi2(iim,jjm+1),airefi2(iim,jjm+1)
+
+          REAL*4 mfu(kon,kev), mfd(kon,kev)
+          REAL*4 en_u(kon,kev), de_u(kon,kev)
+          REAL*4 en_d(kon,kev), de_d(kon,kev)
+          REAL*4 coefh(kon,kev)
+		REAL*4 t(kon,kev)
+
+          REAL*4 mfu2(iim,jjm+1,kev), mfd2(iim,jjm+1,kev)
+          REAL*4 en_u2(iim,jjm+1,kev), de_u2(iim,jjm+1,kev)
+          REAL*4 en_d2(iim,jjm+1,kev), de_d2(iim,jjm+1,kev)
+          REAL*4 coefh2(iim,jjm+1,kev)
+		REAL*4 t2(iim,jjm+1,kev)
+
+          REAL*4 pl(kev)
+          integer irec
+          integer*4 xid,yid,zid,tid
+          integer*4 zrec,zim,zjm
+          integer*4 ncrec,nckon,nckev,ncim,ncjm
+
+          real*4 airefi(kon)
+          character namedim
+
+c  !! attention !!
+c attention il y a aussi le pb de def kon
+c dim de phis??
+
+          REAL*4 frac_impa(kon,kev), frac_nucl(kon,kev)
+          REAL*4 frac_impa2(iim,jjm+1,kev), 
+     .     frac_nucl2(iim,jjm+1,kev)
+          REAL*4 pyu1(kon), pyv1(kon)
+          REAL*4 pyu12(iim,jjm+1), pyv12(iim,jjm+1)
+          REAL*4 ftsol(kon,nbsrf)
+          REAL*4 psrf(kon,nbsrf) 
+	  REAL*4 ftsol1(kon),ftsol2(kon),ftsol3(kon),ftsol4(kon)
+          REAL*4 psrf1(kon),psrf2(kon),psrf3(kon),psrf4(kon)
+          REAL*4 ftsol12(iim,jjm+1),ftsol22(iim,jjm+1),
+     .     ftsol32(iim,jjm+1),
+     .     ftsol42(iim,jjm+1)
+          REAL*4 psrf12(iim,jjm+1),psrf22(iim,jjm+1),psrf32(iim,jjm+1),
+     .     psrf42(iim,jjm+1)
+	
+	  integer ncidp
+          save ncidp
+          integer varidmfu, varidmfd, varidps, varidenu, variddeu	
+		integer varidt
+          integer varidend,varidded,varidch,varidfi,varidfn
+          integer varidyu1,varidyv1,varidpl,varidai,varididvt
+          integer varidfts1,varidfts2,varidfts3,varidfts4
+          integer varidpsr1,varidpsr2,varidpsr3,varidpsr4
+          save varidmfu, varidmfd, varidps, varidenu, variddeu
+		save varidt
+          save varidend,varidded,varidch,varidfi,varidfn
+          save varidyu1,varidyv1,varidpl,varidai,varididvt
+          save varidfts1,varidfts2,varidfts3,varidfts4
+          save varidpsr1,varidpsr2,varidpsr3,varidpsr4
+
+          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
+
+            varidmfu=NCVID(ncidp,'mfu',rcode)
+            print*,'ncidp,varidmfu',ncidp,varidmfu
+
+		varidt=NCVID(ncidp,'t',rcode)
+                print*,'ncidp,varidt',ncidp,varidt
+
+            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
+	
+	    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,nckev)
+            status = nf_inq_dim(ncidp,tid,namedim,ncrec)
+	
+            zrec=ncrec
+            zkev=nckev
+            zim=ncim
+            zjm=ncjm
+	
+	    zkon=zim*(zjm-2)+2
+	
+	    write(*,*) 'read_pstoke : zrec = ', zrec
+            write(*,*) 'read_pstoke : kev = ', zkev
+            write(*,*) 'read_pstoke : zim = ', zim 
+            write(*,*) 'read_pstoke : zjm = ', zjm
+            write(*,*) 'read_pstoke : kon = ', zkon
+
+c niveaux de pression
+
+            status=NF_GET_VARA_REAL(ncidp,varidpl,1,kev,pl)
+
+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 
+c phis
+      status=NF_GET_VARA_REAL(ncidp,varidps,start,count,phisfi2)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,phisfi2,phisfi)
+
+c aire
+      status=NF_GET_VARA_REAL(ncidp,varidai,start,count,airefi2)
+      call gr_ecrit_fi(1,kon,iim,jjm+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)=kev
+      count(4)=1
+
+c frac_impa 
+
+      status=NF_GET_VARA_REAL(ncidp,varidfi,start,count,frac_impa2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,frac_impa2,frac_impa)
+
+c frac_nucl 
+
+      status=NF_GET_VARA_REAL(ncidp,varidfn,start,count,frac_nucl2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,frac_nucl2,frac_nucl)
+
+c abder t
+      status=NF_GET_VARA_REAL(ncidp,varidt,start,count,t2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,t2,t)
+
+c mfu
+      status=NF_GET_VARA_REAL(ncidp,varidmfu,start,count,mfu2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,mfu2,mfu)
+
+c mfd
+      status=NF_GET_VARA_REAL(ncidp,varidmfd,start,count,mfd2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,mfd2,mfd)
+
+c en_u 
+      status=NF_GET_VARA_REAL(ncidp,varidenu,start,count,en_u2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,en_u2,en_u)
+
+c de_u 
+      status=NF_GET_VARA_REAL(ncidp,variddeu,start,count,de_u2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,de_u2,de_u)
+
+c en_d 
+      status=NF_GET_VARA_REAL(ncidp,varidend,start,count,en_d2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,en_d2,en_d)
+
+c de_d 
+      status=NF_GET_VARA_REAL(ncidp,varidded,start,count,de_d2)
+      call gr_ecrit_fi(kev,kon,iim,jjm+1,de_d2,de_d)
+
+c coefh 
+	print*,'LECTURE de coefh a irec =',irec
+       status=NF_GET_VARA_REAL(ncidp,varidch,start,count,coefh2)
+       call gr_ecrit_fi(kev,kon,iim,jjm+1,coefh2,coefh)
+
+      start(3)=irec
+      start(4)=0
+      count(3)=1
+      count(4)=0
+
+c pyu1
+	print*,'LECTURE de yu1 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidyu1,start,count,pyu12)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,pyu12,pyu1)
+
+c pyv1
+        print*,'LECTURE de yv1 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidyv1,start,count,pyv12)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,pyv12,pyv1)
+
+c ftsol1
+        print*,'LECTURE de ftsol1 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidfts1,start,count,ftsol12)
+       call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol12,ftsol1)
+
+c ftsol2
+        print*,'LECTURE de ftsol2 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidfts2,start,count,ftsol22)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol22,ftsol2)
+
+c ftsol3
+	 print*,'LECTURE de ftsol3 a irec =',irec
+      status=NF_GET_VARA_REAL(ncidp,varidfts3,start,count,ftsol32)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol32,ftsol3)
+
+c ftsol4
+      status=NF_GET_VARA_REAL(ncidp,varidfts4,start,count,ftsol42)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,ftsol42,ftsol4)
+
+c psrf1 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr1,start,count,psrf12)
+c      call dump2d(iip1-1,jjm+1,psrf12,'PSRF1NC')
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf12,psrf1)
+
+c psrf2 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr2,start,count,psrf22)
+c      call dump2d(iip1-1,jjm+1,psrf22,'PSRF2NC')
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf22,psrf2)
+
+c psrf3 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr3,start,count,psrf32)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf32,psrf3)
+
+c psrf4 
+      status=NF_GET_VARA_REAL(ncidp,varidpsr4,start,count,psrf42)
+      call gr_ecrit_fi(1,kon,iim,jjm+1,psrf42,psrf4)
+	
+          do i = 1,kon
+	
+	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
+
