program tripgr4D

c--------------------------------------------------------------------

c ***IN***

c But : Lecture d'un fichier grille 4D
c       du champ de T (cree par lectgr)
c	ou de son anomalie par rap. a annee-type (cree par anomalie)
c
c
c	Ttes les grilles en entree sont de la forme:

c        1979-1996 -> it=1,108 (par bimois) 
c	 120E,290E -> ilon=1,35 (pas de 5 deg.)
c	  30S, 30N -> ilat=1,61 (pas de 1 deg.) 
c        0-450m    -> iz=1,34 
c	(pas de 10m de 0 a 250m et 25m de 275m a 450m)

c ***OUT***

c    Creation: 1/ fichier colonne x-plettes:
c 			Long, Lat, val, ect de val 
c		(moy en tps) pour plot avec pplus (csssxy.ppc)
c	       2/ fichier colonne (data_ft) (??)
c			xnumois,prof,valeur,an,mois
c			pour plot avec ft.ppc ou con_tz.ppc
c    	       3/ fichier colonne x-plettes:
c 			Long, Numero de mois, val, ect de val
c		(moy en lat) pour plot avec pplus (csssLt.ppc)
c    	       4/ fichier colonne x-plettes:
c 			Lat, Numero de mois, val, ect de val
c		(moy en long) pour plot avec pplus  cssslt.ppc)
c              5/ fichier colonne x-plettes:
c			Long, Prof, val, ect de val
c		(moy en tps) pour plot avec pplus (nom prog?)
c	       6/ fichier colonne x-plettes:
c			Lat, Prof, val, ect de val
c		(moy en tps) pour plot avec pplus  nom prog?)
c
c       pour des donnees choisies :
c	- entre 2 indices de longitude
c	latitude et mois et pour 1 indice de profondeur (champs
c	horizontaux, longitude-temps ou latitude-temps) 
c	- entre 2 indices de longitude et profondeur et pour 1
c	indice de latitude (champs verticaux zonaux)
c	-entre 2 indices de latitude et profondeur et pour 1 ind.
c	de longitude (champs verticaux meridiens)
c   Compilation : ft tripgr4D
c------------------------------------------------------------------

	parameter(NLO=35,NLA=61,NZ=34,NTI=108)
	dimension rniv(NZ)	
        real val(NLO,NLA,NZ,NTI)
	double precision fsst(NTI)
     	double precision dmiss,valmoy,valect
	character*256 dirin,dirout,namin,namout
    	character*1 iat,ichoix
	integer lnblnk

	data amiss/1.e35/
	dmiss=dble(amiss)

c---------------------------
c	Initialisation
c---------------------------

	do i=1,NLO
	  do j=1,NLA
	    do k=1,NZ
		do it=1,NTI
		  val(i,j,k,it)=amiss
		enddo
	    enddo
	  enddo
	enddo

      rlo1=120.
      rlo2=290.
      rla1=-30.
      rla2= 30.
      palon = (rlo2-rlo1)/float(NLO-1)
      palat = (rla2-rla1)/float(NLA-1)
      write(6,*) 'palon,palat = ',palon,palat

c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c		lecture du fichier AN
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

      dirin='/data/services/archives/temperature/xbt_pacifique/data/'
      write(0,'(//,5x,"Nom grille 4D (GT4D, GT4Derr) : ",$ )')
      read(5,*) namin
      namin=dirin(1:lnblnk(dirin))//namin
      write(6,*) namin
      open (unit=33,file=namin,form='unformatted',status='old')
      read(33) nxi,nyi,nzi,ntime
      read(33) xlon1i,xlon2i,xlat1i,xlat2i,iz1i,iz2i,it1i,it2i
      read(33) (rniv(k),k=iz1i,iz2i)
      write(6,*) nxi,nyi,nzi,ntime
      write(6,*) xlon1i,xlon2i,xlat1i,xlat2i,iz1i,iz2i,it1i,it2i
      write(6,*) (rniv(k),k=iz1i,iz2i)
       
      ilon1i=1+int((xlon1i-rlo1)/palon)
      ilon2i=ilon1i+nxi-1
      ilat1i=1+int((xlat1i-rla1)/palat)
      ilat2i=ilat1i+nyi-1
      write(6,*) "ilon1i,ilon2i,ilat1i,ilat2i =", 
     &            ilon1i,ilon2i,ilat1i,ilat2i

      read(33) ((((val(i,j,k,it),i=ilon1i,ilon2i),j=ilat1i,ilat2i),
     & k=iz1i,iz2i),it=it1i,it2i)
      close(33)

c---------------------------------------
c	Passe bande (pour val. absentes)
c---------------------------------------

	tmin=-20.
	tmax= 30.

	do i=1,NLO
	  do j=1,NLA
	    do k=1,NZ
		do it=1,NTI
		  if(val(i,j,k,it).ge.tmax.or.val(i,j,k,it).le.tmin) then
		    val(i,j,k,it)=amiss
		  endif
		enddo
	    enddo
	  enddo
	enddo


c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      write(0,*) ''
      write(0,*) '  Choix des indices de longitude, latitude, prof '
      write(0,*) ' et temps ou sont calcules moyennes et ecart-types'
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

88    continue

      write(0,'(//"*** SELECTION LATITUDES ***")')
        write(0,'(10x"Latitudes MIN, MAX (-30,30,1): ",$)')
        read(5,*) xlat1,xlat2
        write(0,*) xlat1,xlat2
        ilat1=1+int((xlat1-rla1)/palat)        
        ilat2=1+int((xlat2-rla1)/palat)
	write(0,*) ilat1,ilat2,' = ilat1,ilat2'

      write(0,'(//"*** SELECTION LONGITUDES ***")')
        write(0,'(10x"Longitudes MIN, MAX (120,290,5): ",$)')
        read(5,*) xlon1,xlon2
	if(xlon1.lt.0.) xlon1=xlon1+360.
	if(xlon2.lt.0.) xlon2=xlon2+360.
        write(0,*) xlon1,xlon2
        ilon1=1+int((xlon1-rlo1)/palon)        
        ilon2=1+int((xlon2-rlo1)/palon)
	write(0,*) ilon1,ilon2,' = ilon1,ilon2'

      write(0,'(//"*** SELECTION PROFONDEUR ***")')
      
        write(0,'(10x,"Affichage Indice profondeur, o/n : ",$)')
	read(5,*) iat
	if(iat.eq.'o') then
	  do iz=iz1i,iz2i
	    write(0,'(i2,2x,f5.1)') iz,rniv(iz)
	  enddo
	endif

        write(0,'(/,10x," Indices profondeur Debut,Fin : ",$)')
        read(5,*) iz1,iz2
	nbz=iz2-iz1+1
	write(0,*) nbz,' = nbz'


      write(0,'(//"*** SELECTION TEMPS (mois) ***")')
        write(0,'(10x,"Affichage ecran des numeros bimois o/n : ",$)')
	read(5,*) iat
	if(iat.eq.'o') then
	  ia1=79
	  ia2=96
	  do ia=ia1,ia2
	    write(0,'(i2,2x,12i4)') ia,((ia-ia1)*6+it,it=1,6)
	  enddo
	endif
	  
        write(0,'(/,10x," No bimois Debut, Fin : ",$)')
        read(5,*) it1,it2

  


c--------------------------------------------------------------
      write(0,'(//"*** NATURE FICHIER OUT ***")')

        write(0,'(/5x,"OUT, Long, Lat, Val, Ect,")')
        write(0,'(5x,"   et fichier climatologique    entrer 1 : ")')
        write(0,'(5x,"OUT, Xmois,Prof,Valeur,An,Bimois  entrer 2 : ")')
        write(0,'(5x,"OUT, Long, temps, Valeur,       entrer 3 : ")')
        write(0,'(5x,"OUT, Lat, temps, Valeur,        entrer 4 : ")')
        write(0,'(5x,"OUT, Long, prof, Valeur,        entrer 5 : ")')
        write(0,'(5x,"OUT, Lat, prof, Valeur,         entrer 6 : ",$)')
        read(5,*) iout
c--------------------------------------------------------------

      write(0,'(/30x,"Facteur multiplicatif sur Val, Ect: ",$)')
      read(5,*) xmul
      dirout='./'
c
c-> IOUT=1

      if(iout.eq.1) then
        write(0,'(//,5x," Fichier X-plettes en Sortie : ",$)')
        read(*,'(a)') namout
	namout=dirout(1:lnblnk(dirout))//namout
        open(unit=34,file=namout,form='formatted',status='unknown')

        do i=ilon1,ilon2
	   xlon = rlo1 + float(i-1)*palon 
           do j=ilat1,ilat2
	      xlat = rla1 + float(j-1)*palat
	      itt=0 
              do it=it1,it2
		  itt=itt+1
		  fsst(itt)=dble(val(i,j,iz1,it))
	      enddo
	      call dmoyect(fsst,itt,dmiss,valmoy,valect)
	      if(valmoy.lt.dmiss) then
		xvalmoy=sngl(valmoy)*xmul
	        xect=sngl(valect)*xmul
	      else
		xvalmoy=sngl(valmoy)
	        xect=sngl(valect)
	      endif
	      write(34,77) xlon,xlat,xvadmoy,xect
c	      write(0,*) xlon,xlat,xvalmoy,xect

	   enddo
        enddo

        close(34)

      endif 
c
c-> IOUT=2

      if(iout.eq.2) then
        write(0,'(//,5x," Fichier X-plettes en Sortie : ",$)')
        read(*,'(a)') namout
	namout=dirout(1:lnblnk(dirout))//namout
        open(unit=34,file=namout,form='formatted',status='unknown')

	do it=it1,it2
	 do iz=iz1,iz2

	  xx=float(it)/6.
	  ix=int(it/6.)
	  if(xx.eq.float(ix)) ix=ix-1
	  ian=1979+ix
	  ibimois=mod(it,6)
	  if(ibimois.eq.0) ibimois=6
c	  write(6,*) it,ian,ibimois

	  xsom=0.
	  xmoy=0.
	  xnobs=0.
          do i=ilon1,ilon2
             do j=ilat1,ilat2
	       if(val(i,j,iz,it).ne.amiss) then
		 xsom=xsom+val(i,j,iz,it)
		 xnobs=xnobs+1.
	       endif
	     enddo
	  enddo

	  if(xnobs.gt.0.) then
	    xmoy=xsom/xnobs
	  else
	    xmoy=amiss
	  endif

	  if(xmoy.lt.amiss) then
	    xmoy2=xmoy*xmul
	  else
	    xmoy2=xmoy
	  endif

	  write(34,77) float(it),rniv(iz),xmoy2,float(ian),float(ibimois)
c	  write(0,*) float(it),rniv(iz),xmoy2,float(ian),float(imois)

 	 enddo
	enddo

        close(34)
      endif
c
c-> IOUT=3

      if(iout.eq.3) then
        write(0,'(//,5x," Fichier X-plettes en Sortie : ",$)')
        read(*,'(a)') namout
	namout=dirout(1:lnblnk(dirout))//namout
        open(unit=34,file=namout,form='formatted',status='unknown')
c
        do i=ilon1,ilon2
	  xlon = rlo1 + float(i-1)*palon 
            do it=it1,it2
	      ila=0 
              do j=ilat1,ilat2
		  ila=ila+1
		  fsst(ila)=dble(val(i,j,iz1,it))
	      enddo
	      call dmoyect(fsst,ila,dmiss,valmoy,valect)
	      if(valmoy.lt.dmiss) then
		xvalmoy=sngl(valmoy)*xmul
	        xect=sngl(valect)*xmul
	      else
		xvalmoy=sngl(valmoy)
	        xect=sngl(valect)
	      endif
	      write(34,77) xlon,float(it),xvalmoy,xect
c	      write(0,*) xlon,float(it),xvalmoy,xect
	   enddo
        enddo
        close(34)
      endif 
c
c-> IOUT=4

      if(iout.eq.4) then
        write(0,'(//,5x," Fichier X-plettes en Sortie : ",$)')
        read(*,'(a)') namout
	namout=dirout(1:lnblnk(dirout))//namout
        open(unit=34,file=namout,form='formatted',status='unknown')
c
         do j=ilat1,ilat2
	    xlat = rla1 + float(j-1)*palat
            do it=it1,it2
	      ilo=0 
              do i=ilon1,ilon2
		  ilo=ilo+1
		  fsst(ilo)=dble(val(i,j,iz1,it))
	      enddo
	      call dmoyect(fsst,ilo,dmiss,valmoy,valect)
	      if(valmoy.lt.dmiss) then
		xvalmoy=sngl(valmoy)*xmul
	        xect=sngl(valect)*xmul
	      else
		xvalmoy=sngl(valmoy)
	        xect=sngl(valect)
	      endif
	      write(34,77) xlat,float(it),xvalmoy,xect
c	      write(0,*) xlat,float(it),xvalmoy,xect
	   enddo
        enddo
        close(34)
      endif 
c
c-> IOUT=5

      if(iout.eq.5) then
        write(0,'(//,5x," Fichier X-plettes en Sortie : ",$)')
        read(*,'(a)') namout
	namout=dirout(1:lnblnk(dirout))//namout
        open(unit=34,file=namout,form='formatted',status='unknown')

        do i=ilon1,ilon2
	   xlon = rlo1 + float(i-1)*palon 
           do k=iz1,iz2
	      xprof=rniv(k)
	      itt=0 
              do it=it1,it2
		  itt=itt+1
		  fsst(itt)=dble(val(i,ilat1,k,it))
	      enddo
	      call dmoyect(fsst,itt,dmiss,valmoy,valect)
	      if(valmoy.lt.dmiss) then
		xvalmoy=sngl(valmoy)*xmul
	        xect=sngl(valect)*xmul
	      else
		xvalmoy=sngl(valmoy)
	        xect=sngl(valect)
	      endif
	      write(34,77) xlon,xprof,xvalmoy,xect
c	      write(0,*) xlon,xprof,xvalmoy,xect
	   enddo
        enddo

        close(34)

      endif 

c
c-> IOUT=6

      if(iout.eq.6) then
        write(0,'(//,5x," Fichier X-plettes en Sortie : ",$)')
        read(*,'(a)') namout
	namout=dirout(1:lnblnk(dirout))//namout
        open(unit=34,file=namout,form='formatted',status='unknown')

        do j=ilat1,ilat2
	   xlat = rla1 + float(j-1)*palat 
           do k=iz1,iz2
	      xprof=rniv(k)
	      itt=0 
              do it=it1,it2
		  itt=itt+1
		  fsst(itt)=dble(val(ilon1,j,k,it))
	      enddo
	      call dmoyect(fsst,itt,dmiss,valmoy,valect)
	      if(valmoy.lt.dmiss) then
		xvalmoy=sngl(valmoy)*xmul
	        xect=sngl(valect)*xmul
	      else
		xvalmoy=sngl(valmoy)
	        xect=sngl(valect)
	      endif
	      write(34,77) xlat,xprof,xvalmoy,xect
c	      write(0,*) xlat,xprof,xvalmoy,xect
	   enddo
        enddo

        close(34)

      endif 


      write(0,'(//,5x," Autres Choix o/n : ",$)')
      read(5,'(a)') ichoix
      if(ichoix.eq.'o') go to 88


999     stop 000
77	format(6e13.4)
        end

c---------------------------------------------------
        subroutine dmoyect(dx,np,dmiss,dxmoy,dxrms)
c        
c       A partir d'un tableau dx(np) ou les valeurs
c       absentes st a "dmiss", calcul moyennne et
c       ecart-type
c--------------------------------------------------

        double precision dx(np),dsx,dsx2,dxmoy,dxrms,dmiss

        nvalid=0
        dsx=0.d0
        dsx2=0.d0
        do i=1,np
          if(dx(i).ne.dmiss) then
            dsx = dsx + dx(i)
            dsx2 = dsx2 + dx(i)*dx(i)
            nvalid = nvalid+1
          endif
        enddo
         
        if(nvalid.ge.1) then
          dxmoy=dsx/nvalid
	  if(nvalid.eq.1) then
	    dxrms=dmiss
	  else
            dxrms=(nvalid*dsx2-dsx*dsx) / (nvalid*(nvalid-1))
            dxrms=dsqrt(dxrms)
	  endif
        else
          dxmoy=5*dmiss
          dxrms=5*dmiss
        endif

        return
        end