!
! $Id: grid_noro.F 1999 2014-03-20 09:57:19Z aborella $
!
c
c
      SUBROUTINE grid_noro(imdep, jmdep, xdata, ydata, zdata,
     .             imar, jmar, x, y,
     .             zphi,zmea,zstd,zsig,zgam,zthe,
     .             zpic,zval,mask)
c=======================================================================
c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
c
c      Compute the Parameters of the SSO scheme as described in
c      LOTT & MILLER (1997) and LOTT(1999).
c      Target points are on a rectangular grid:
c      iim+1 latitudes including North and South Poles;
c      jjm+1 longitudes, with periodicity jjm+1=1.
c      aux poles.  At the poles the fields value is repeated
c      jjm+1 time.
c      The parameters a,b,c,d represent the limite of the target
c      gridpoint region. The means over this region are calculated
c      from USN data, ponderated by a weight proportional to the 
c      surface occupated by the data inside the model gridpoint area.
c      In most circumstances, this weight is the ratio between the
c      surface of the USN gridpoint area and the surface of the
c      model gridpoint area. 
c
c           (c)
c        ----d-----
c        | . . . .|
c        |        |
c     (b)a . * . .b(a)
c        |        |
c        | . . . .|
c        ----c-----
c           (d)
C=======================================================================
c INPUT:
c        imdep, jmdep: dimensions X and Y input field
c        xdata, ydata: coordinates X and Y input field
c        zdata: Input field
c        In this version it is assumed that the entry data come from
c        the USNavy dataset: imdep=iusn=2160, jmdep=jusn=1080.
c OUTPUT:
c        imar, jmar: dimensions X and Y Output field
c        x, y: ccordinates  X and Y Output field.
c             zmea:  Mean orographie   
c             zstd:  Standard deviation
c             zsig:  Slope
c             zgam:  Anisotropy
c             zthe:  Orientation of the small axis
c             zpic:  Maximum altitude
c             zval:  Minimum altitude
C=======================================================================

      IMPLICIT INTEGER (I,J)
      IMPLICIT REAL(X,Z) 
      
	  parameter(iusn=2160,jusn=1080,iext=216, epsfra = 1.e-5)
#include "dimensions.h"
	  REAL xusn(iusn+2*iext),yusn(jusn+2)	
      REAL zusn(iusn+2*iext,jusn+2)

      INTEGER imdep, jmdep
      REAL xdata(imdep),ydata(jmdep) 
      REAL zdata(imdep,jmdep)
c
      INTEGER imar, jmar
  
C INTERMEDIATE FIELDS  (CORRELATIONS OF OROGRAPHY GRADIENT)

      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
      REAL weight(iim+1,jjm+1)

C CORRELATIONS OF USN OROGRAPHY GRADIENTS

      REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2)
      REAL zxtzyusn(iusn+2*iext,jusn+2)
      REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
      REAL zmea0(imar+1,jmar) ! GK211005 (CG)
      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
      REAL zpic(imar+1,jmar),zval(imar+1,jmar)
cxxx PB     integer mask(imar+1,jmar)
      real mask(imar+1,jmar), mask_tmp(imar+1,jmar)
      real num_tot(2200,1100),num_lan(2200,1100)
c
      REAL a(2200),b(2200),c(1100),d(1100)
      logical masque_lu
c
      print *,' parametres de l orographie a l echelle sous maille' 
      xpi=acos(-1.)
      rad    = 6 371 229.
      zdeltay=2.*xpi/REAL(jusn)*rad
c
c utilise-t'on un masque lu?
c
      masque_lu = .true.
      if (maxval(mask) == -99999 .and. minval(mask) == -99999) then
        masque_lu= .false.
        masque = 0.0
      endif
      write(*,*)'Masque lu', masque_lu
c
c  quelques tests de dimensions:
c    
c
      if(iim.ne.imar) STOP 'Problem dim. x'
      if(jjm.ne.jmar-1) STOP 'Problem dim. y'
      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
         PRINT*, 'imar or jmar too big', imar, jmar
         CALL ABORT_GCM("GRID_NORO", "", 1)
      ENDIF

      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
         print *,' imdep or jmdep bad dimensions:',imdep,jmdep
         call abort_gcm("grid_noro", "", 1)
      ENDIF

      IF(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
        print *,' imar or jmar bad dimensions:',imar,jmar
        call abort_gcm("grid_noro", "", 1)
      ENDIF


c      print *,'xdata:',xdata
c      print *,'ydata:',ydata
c      print *,'x:',x
c      print *,'y:',y
c
C  EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
C  BOUNDARIES:
c
      DO j=1,jusn
        yusn(j+1)=ydata(j)
      DO i=1,iusn
        zusn(i+iext,j+1)=zdata(i,j)
        xusn(i+iext)=xdata(i)
      ENDDO
      DO i=1,iext
        zusn(i,j+1)=zdata(iusn-iext+i,j)
        xusn(i)=xdata(iusn-iext+i)-2.*xpi
        zusn(iusn+iext+i,j+1)=zdata(i,j)
        xusn(iusn+iext+i)=xdata(i)+2.*xpi
      ENDDO
      ENDDO

        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
        yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
       DO i=1,iusn/2+iext
        zusn(i,1)=zusn(i+iusn/2,2)
        zusn(i+iusn/2+iext,1)=zusn(i,2)
        zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
        zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
       ENDDO
c  
c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
C     ( REGULAR GRID)
c
      a(1) = x(1) - (x(2)-x(1))/2.0
      b(1) = (x(1)+x(2))/2.0
      DO i = 2, imar
         a(i) = b(i-1)
         b(i) = (x(i)+x(i+1))/2.0
      ENDDO
      a(imar+1) = b(imar)
      b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/2.0

      c(1) = y(1) - (y(2)-y(1))/2.0
      d(1) = (y(1)+y(2))/2.0
      DO j = 2, jmar-1
         c(j) = d(j-1)
         d(j) = (y(j)+y(j+1))/2.0
      ENDDO
      c(jmar) = d(jmar-1)
      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
c
c  initialisations:
c
      DO i = 1, imar+1
      DO j = 1, jmar
         weight(i,j) = 0.0
         zxtzx(i,j)  = 0.0
         zytzy(i,j)  = 0.0
         zxtzy(i,j)  = 0.0
         ztz(i,j)    = 0.0
         zmea(i,j)   = 0.0
         zpic(i,j)  =-1.E+10
         zval(i,j)  = 1.E+10
      ENDDO
      ENDDO
c
c  COMPUTE SLOPES CORRELATIONS ON USN GRID
c
         DO j = 1,jusn+2 
         DO i = 1, iusn+2*iext
            zytzyusn(i,j)=0.0
            zxtzxusn(i,j)=0.0
            zxtzyusn(i,j)=0.0
         ENDDO
         ENDDO


         DO j = 2,jusn+1 
            zdeltax=zdeltay*cos(yusn(j))
         DO i = 2, iusn+2*iext-1
            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
         ENDDO
         ENDDO
c
c  SUMMATION OVER GRIDPOINT AREA
c 
      zleny=xpi/REAL(jusn)*rad
      xincr=xpi/2./REAL(jusn)
       DO ii = 1, imar+1
       DO jj = 1, jmar
       num_tot(ii,jj)=0.
       num_lan(ii,jj)=0.
c        PRINT *,' iteration ii jj:',ii,jj
         DO j = 2,jusn+1 
c         DO j = 3,jusn 
            zlenx=zleny*cos(yusn(j))
            zdeltax=zdeltay*cos(yusn(j))
            zbordnor=(c(jj)-yusn(j)+xincr)*rad
            zbordsud=(yusn(j)-d(jj)+xincr)*rad
            weighy=AMAX1(0.,
     *             amin1(zbordnor,zbordsud,zleny))
         IF(weighy.ne.0)THEN
         DO i = 2, iusn+2*iext-1
            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
            weighx=AMAX1(0.,
     *             amin1(zbordest,zbordoue,zlenx))
            IF(weighx.ne.0)THEN
            num_tot(ii,jj)=num_tot(ii,jj)+1.0
            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
            weight(ii,jj)=weight(ii,jj)+weighx*weighy
            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
c mean
            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
c peacks
            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
c valleys
            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
            ENDIF
         ENDDO
         ENDIF
         ENDDO
       ENDDO
       ENDDO
c
c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
C  LOTT (1999) SSO SCHEME.
c
      zllmmea=0.
      zllmstd=0.
      zllmsig=0.
      zllmgam=0.
      zllmpic=0.
      zllmval=0.
      zllmthe=0.
      zminthe=0.
c     print 100,' '
c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH') 
       DO ii = 1, imar+1
       DO jj = 1, jmar
         IF (weight(ii,jj) .NE. 0.0) THEN
c  Mask
cXXX           if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
cXXX             mask(ii,jj)=1
cXXX           else
cXXX             mask(ii,jj)=0
cXXX           ENDIF
             if (.not. masque_lu) then
               mask(ii,jj) = num_lan(ii,jj)/num_tot(ii,jj)
             endif
c  Mean Orography:
           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
c  Standard deviation:
           zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
         ELSE
            PRINT*, 'probleme,ii,jj=', ii,jj
         ENDIF
       ENDDO
       ENDDO

C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:

       DO ii = 1, imar+1
         zxtzx(ii,1)=zxtzx(ii,2)
         zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
         zxtzy(ii,1)=zxtzy(ii,2)
         zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
         zytzy(ii,1)=zytzy(ii,2)
         zytzy(ii,jmar)=zytzy(ii,jmar-1)
       ENDDO

C  FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.

C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.

       zmea0(:,:) = zmea(:,:) ! GK211005 (CG) on sauvegarde la topo non lissee
       CALL MVA9(zmea,iim+1,jjm+1)
       CALL MVA9(zstd,iim+1,jjm+1)
       CALL MVA9(zpic,iim+1,jjm+1)
       CALL MVA9(zval,iim+1,jjm+1)
       CALL MVA9(zxtzx,iim+1,jjm+1)
       CALL MVA9(zxtzy,iim+1,jjm+1) 
       CALL MVA9(zytzy,iim+1,jjm+1)
CXXX   Masque prenant en compte maximum de terre
CXXX  On seuil a 10% de terre de terre car en dessous les parametres de surface n'on
CXXX pas de sens (PB)
       mask_tmp= 0.0
       WHERE(mask .GE. 0.1) mask_tmp = 1.

       DO ii = 1, imar
       DO jj = 1, jmar
         IF (weight(ii,jj) .NE. 0.0) THEN
c  Coefficients K, L et M:
           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
           xm=zxtzy(ii,jj)
           xp=xk-sqrt(xl**2+xm**2)
           xq=xk+sqrt(xl**2+xm**2)
           xw=1.e-8
           if(xp.le.xw) xp=0.
           if(xq.le.xw) xq=xw
           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
c slope: 
cXXX           zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
cXXXc isotropy:
cXXX           zgam(ii,jj)=xp/xq*mask(ii,jj)
cXXXc angle theta:
cXXX           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask(ii,jj)
cXXX           zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
cXXX           zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
cXXX           zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
cXXX           zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
cXXX           zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
CXX* PB modif pour maque de terre fractionnaire
c slope: 
           zsig(ii,jj)=sqrt(xq)*mask_tmp(ii,jj)
c isotropy:
           zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
c angle theta:
           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask_tmp(ii,jj)
           ! GK211005 (CG) ne pas forcement lisser la topo
           ! zphi(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
           zphi(ii,jj)=zmea0(ii,jj)*mask_tmp(ii,jj)
           !
           zmea(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
           zpic(ii,jj)=zpic(ii,jj)*mask_tmp(ii,jj)
           zval(ii,jj)=zval(ii,jj)*mask_tmp(ii,jj)
           zstd(ii,jj)=zstd(ii,jj)*mask_tmp(ii,jj)
c          print 101,ii,jj,
c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
c    *           zthe(ii,jj)
c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
         ELSE
c           PRINT*, 'probleme,ii,jj=', ii,jj
         ENDIF
      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
      zminthe=amin1(zthe(ii,jj),zminthe)
      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
      zllmval=AMAX1(zval(ii,jj),zllmval)
       ENDDO
       ENDDO
      print *,'  MEAN ORO:',zllmmea
      print *,'  ST. DEV.:',zllmstd
      print *,'  PENTE:',zllmsig
      print *,' ANISOTROP:',zllmgam
      print *,'  ANGLE:',zminthe,zllmthe	
      print *,'  pic:',zllmpic
      print *,'  val:',zllmval
      
C
c gamma and theta a 1. and 0. at poles
c
      DO jj=1,jmar
      zmea(imar+1,jj)=zmea(1,jj)
      zphi(imar+1,jj)=zphi(1,jj)
      zpic(imar+1,jj)=zpic(1,jj)
      zval(imar+1,jj)=zval(1,jj)
      zstd(imar+1,jj)=zstd(1,jj)
      zsig(imar+1,jj)=zsig(1,jj)
      zgam(imar+1,jj)=zgam(1,jj)
      zthe(imar+1,jj)=zthe(1,jj)
      ENDDO


      zmeanor=0.0
      zmeasud=0.0
      zstdnor=0.0
      zstdsud=0.0
      zsignor=0.0
      zsigsud=0.0
      zweinor=0.0
      zweisud=0.0
      zpicnor=0.0
      zpicsud=0.0                                   
      zvalnor=0.0
      zvalsud=0.0 

      DO ii=1,imar
      zweinor=zweinor+              weight(ii,   1)
      zweisud=zweisud+              weight(ii,jmar)
      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
      ENDDO

      DO ii=1,imar+1
      zmea(ii,   1)=zmeanor/zweinor
      zmea(ii,jmar)=zmeasud/zweisud
      zphi(ii,   1)=zmeanor/zweinor
      zphi(ii,jmar)=zmeasud/zweisud
      zpic(ii,   1)=zpicnor/zweinor
      zpic(ii,jmar)=zpicsud/zweisud
      zval(ii,   1)=zvalnor/zweinor
      zval(ii,jmar)=zvalsud/zweisud
      zstd(ii,   1)=zstdnor/zweinor
      zstd(ii,jmar)=zstdsud/zweisud
      zsig(ii,   1)=zsignor/zweinor
      zsig(ii,jmar)=zsigsud/zweisud
      zgam(ii,   1)=1.
      zgam(ii,jmar)=1.
      zthe(ii,   1)=0.
      zthe(ii,jmar)=0.
      ENDDO

      RETURN
      END

      SUBROUTINE MVA9(X,IMAR,JMAR)

C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS

      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
      real WEIGHTpb(-1:1,-1:1)


      SUM=0.
      DO IS=-1,1
        DO JS=-1,1
          WEIGHTpb(IS,JS)=1./REAL((1+IS**2)*(1+JS**2))
          SUM=SUM+WEIGHTpb(IS,JS)
        ENDDO
      ENDDO
      
c     WRITE(*,*) 'MVA9 ', IMAR, JMAR
c     WRITE(*,*) 'MVA9 ', WEIGHTpb
c     WRITE(*,*) 'MVA9 SUM ', SUM
      DO IS=-1,1
        DO JS=-1,1
          WEIGHTpb(IS,JS)=WEIGHTpb(IS,JS)/SUM
        ENDDO
      ENDDO

      DO J=2,JMAR-1
        DO I=2,IMAR-1
          XF(I,J)=0.
          DO IS=-1,1
            DO JS=-1,1
              XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHTpb(IS,JS)
            ENDDO
          ENDDO
        ENDDO
      ENDDO

      DO J=2,JMAR-1
        XF(1,J)=0.
        IS=IMAR-1
        DO JS=-1,1 
          XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHTpb(-1,JS)
        ENDDO
        DO IS=0,1 
          DO JS=-1,1 
            XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHTpb(IS,JS)
          ENDDO
        ENDDO
        XF(IMAR,J)=XF(1,J)
      ENDDO

      DO I=1,IMAR
        XF(I,1)=XF(I,2)
        XF(I,JMAR)=XF(I,JMAR-1)
      ENDDO

      DO I=1,IMAR
        DO J=1,JMAR
          X(I,J)=XF(I,J)
        ENDDO
      ENDDO

      RETURN
      END