C------------------------------------------------------------
C
C     CODE TO PRODUCE ETA SOIL TEXTURE;                
C     PROGRAM WRITTEN USING THE STRUCTURE OF mountnew.f
C     WRITTEN BY F. MESINGER, W. COLLINS AND Z. JANJIC.
C     (MODIFIED BY S. NICKOVIC FOR SURFACE TYPES)
C
C************************************************************
C
C     PROGRAM HISTORY LOG:
C      Jan 1996 S. NICKOVIC     - ORIGINATOR
C      Jul 1998 A. PAPADOPOULOS - CORRECTION AT COASTAL AREAS
C      Jun 1999 P. KATSAFADOS   - CROSSING 180W
C************************************************************
C                                                           *
C     DATA SOURCE: UNEP/FAO SOIL TYPES                      *
C     WITH 134 CATHEGORIES; RESOLUTION: 2'X2'DEG            *
C     DATA ARRAY STARTS FROM NBD=-90, EBD=-180              *
C                                                           *
C************************************************************
C                                                           *
C     METHOD: THE LARGEST PERCENT OF A SOIL TYPE  THAT      *
C     ENTERS IN THE H GRID-BOX IS USED AS                   * 
C     THE REPRESENTATIVE ONE                                *    
C                                                           *
C************************************************************
C                                                           *
#include "mountnew.inc"
C---------------------
      INCLUDE "textparm.inc"
C---------------------
      PARAMETER(NTEX=10,DLD=1./30.)
      PARAMETER(ISOISEA=10, PERSEA=0.75)  ! 0.75
      PARAMETER (PI=3.1415926,PI2=2.*PI)
      DIMENSION NDB(IMJM,NTEX), NVP(IMJM), PER(NTEX)
      DIMENSION IST(IROW)

      DIMENSION ISOI(IMJM)
      LOGICAL LAST, INSIDE
      DATA DTR/.01745329/
C
      INCLUDE "mount_nl.nml"
      INCLUDE "soil_nl.nml"
C
      IF (ICROSS.EQ.1) THEN
         WRITE(*,*) 'texteta: Crossing 180deg...'
      ELSE
         WRITE(*,*) 'texteta: Normal soil'
      ENDIF 
C------------------------------------------------------------
      OPEN (UNIT=11,FILE='../../../../namelist/name.list'
     &,     FORM='FORMATTED')
      READ(11,MOUNT_NL)
      READ(11,SOIL_NL)
      CLOSE(11)
C------------------------------------------------------------
C
      write(*,*)'IM=',IM,' JM=',JM
      WRITE(*,*)'JCOL=',JCOL,'IROW=',IROW

      DO K=1,IMJM
      ISOI(K)=0
      NVP(K)=0
      DO I=1,NTEX  
      NDB(K,I)=0
      ENDDO
      ENDDO
C
C     HERE STARTS PROCESSING ROW BY ROW
C
      LAST=.FALSE.
      NRRD=0
C
      OPEN (UNIT=99,FILE=FN_TEXT4ETA,FORM='unformatted')
  102 CONTINUE
C
       DO I=1,IROW
       READ(99)IST(I)
       END DO
C
      ALT0D=SOILNBD-(NRRD-1)*DLD   - DLD/2  
      IF (ALT0D.GT.SOILNBD) ALT0D = SOILNBD  !**Check for the NPole** 
      NRRD=NRRD+1
      ILAST=0
C
      IF (NRRD.EQ.JCOL) LAST=.TRUE.
C
      DO 106 ID=1,IROW
C
C       Soil texture at this point is IST(ID)
C
      KSTXT=IST(ID)
C
c************check for crossing 180deg*******************
      IF (ICROSS.EQ.1) THEN
         ALM=(SOILEBD+(ID-1)*DLD)*DTR
         APH=(SOILNBD-(NRRD-1)*DLD)*DTR
         IF (ALM.GT.PI) ALM = ALM - PI2
      ELSEIF (ICROSS.EQ.0) THEN
         ALM=(SOILWBD+(ID-1)*DLD)*DTR
         APH=(SOILNBD-(NRRD-1)*DLD)*DTR
      ENDIF

      CALL PQK(ALM,APH,INSIDE,P,Q,K)
      IF (.NOT.INSIDE) GO TO 106
      ILAST=ILAST+1
C
C       THE SOILTEXTURE (KSTXT) DATA POINT IS WITHIN THE BOX 
C       FORMED BY FOUR NEIGHBORING H POINTS WITH H(K) POINT 
C       AS SOUTH POINT 
C
C                          H(K+NINC)
C
C                    H(K+KNW)      H(K+KNE)
C
C                            H(K)
C
C       NOW CALCULATE THE NUMBER OF POINTS FOR EACH OF 73 
C       SOITYPES THAT ARE INSIDE THE H GRID-BOX:
C
C
           KK=K
           IF (P.GT.0.5 .AND. Q.LE.0.5) KK=K+KNE
           IF (P.GT.0.5 .AND. Q.GT.0.5) KK=K+NINC
           IF (P.LE.0.5 .AND. Q.GT.0.5) KK=K+KNW

           NDB(KK,KSTXT)=NDB(KK,KSTXT)+1
           NVP(KK) = NVP(KK) + 1
C
  106 CONTINUE
      LAST=LAST.OR.(ILAST.EQ.0.AND.ALT0D.LT.TPH0D)
      IF (.NOT.LAST) GO TO 102
C---------------------------------------------------------
C     ALL ROWS OF SOILTEXTURE DATA HAVE BEEN PROCESSED.
C     NOW FIND SOILTEXTURE WHICH HAS MAX CONTRIBUTION IN A 
C     H-GRID BOX; MULTIPLY ISOI WITH SEA-MASK FACTOR TO BE 
C     SURE THAT SEA POINTS ARE SEEN AS SEA POINTS IN 
C     THE SOI DATA
C---------------------------------------------------------
      DO K=1,IMJM
         INIMAX=0
         IINDEX=0
         INIMAX2=0
         IINDEX2=0
         DO I=1,NTEX
            IF (NDB(K,I).GE.INIMAX) THEN
              INIMAX=NDB(K,I)
              IINDEX=I
            ENDIF
         ENDDO
         ISOI(K)=IINDEX
         DO I=1,NTEX
            PER(I) = FLOAT(NDB(K,I))/FLOAT(NVP(K))
         ENDDO
         IF(PER(ISOISEA).LT.PERSEA .AND. IINDEX.EQ.ISOISEA) THEN
           DO I=1,NTEX
             IF (NDB(K,I) .NE. NDB(K,ISOISEA)) THEN
                IF (NDB(K,I).GE.INIMAX2) THEN
                   INIMAX2=NDB(K,I)
                   IINDEX2=I
                ENDIF
             ENDIF
           ENDDO
           ISOI(K)=IINDEX2
cpk           print*,'IINDEX, IINDEX2 :',IINDEX, IINDEX2,NVP(K)
         ENDIF
      ENDDO 
C
      OPEN(UNIT=23,FILE='../tmp/texteta.dat',STATUS='UNKNOWN')
      DO K=1,IMJM
      WRITE(23,*)K,ISOI(K)
      END DO
      CLOSE (23)
C
      OPEN(UNIT=23,FILE=FN_SOIL
     &        ,STATUS='UNKNOWN',FORM='UNFORMATTED')
      WRITE(23)ISOI
      CLOSE (23)
C
      STOP
      END
C***********************************************************************
                             SUBROUTINE PQK
     &(                    ALM,APH,INSIDE,P,Q,K      )
C     ******************************************************************
C     *                                                                *
C     *  ROUTINE TO TRANSFORM ALAMBDA,APHI (LONGITUDE,LATITUDE) INTO   *
C     *  TRANSFORMED (ROTATED) LONGITUDE,LATITUDE, CHECK WHETHER THE   *
C     *  VALUES OBTAINED ARE WITHIN THE MODEL REGION, AND, IF THEY     *
C     *  ARE, CALCULATE VALUES OF P,Q (COORDINATES WITHIN A SQUARE     *
C     *  FORMED BY CONNECTING FOUR NEIGHBORING HEIGHT POINTS) AND K.   *
C     *                                                                *
C     ******************************************************************
C     *                                                                *
C     *  GRID CONSTANTS:                                               *
C     *  WBD,SBD - WESTERN AND SOUTHERN BOUNDARIES ON LL OR TLL GRID   *
C     *            IN DEGREES                                          *
C     *  TLM0D,TPH0D - ANGLES OF ROTATION OF THE LL COORDINATE SYSTEM  *
C     *            IN THE DIRECTION OF LAMBDA AND PHI RESPECTIVELY     *
C     *            IN ORDER TO OBTAIN THE TLL COORDINATE SYSTEM        *
C     *  DLMD,DPHD - MESH SIDES IN DEGREES                             *
C     *                                                                *
C     ******************************************************************
C
#include "mountnew.inc"
                          L O G I C A L
     &  INSIDE
C
                             D A T A
     &  DTR/.01745329/, D50/0.5/, H1/1./
C
C---------------------------CONSTANTS-----------------------------------
      RIM1=IM1
      TLM0=TLM0D*DTR
      TPH0=TPH0D*DTR
      DLM=DLMD*DTR
      DPH=DPHD*DTR
      RDLM=H1/DLM
      RDPH=H1/DPH
      ALMWB=WBD*DTR !+DLM
      APHSB=SBD*DTR !+DPH
      ALMEB=ALMWB   +2*(IM-1)*DLM
      APHNB=APHSB   +  (JM-1)*DPH
      WB=ALMWB !-DLM
      SB=APHSB !-DPH
          STPH0= SIN(TPH0)
          CTPH0= COS(TPH0)
C
C--------------ALM,APH IS LON,LAT IN RADIAN MEASURE;--------------------
C------------TRANSFORM INTO ROTATED LONGITUDE,LATITUDE------------------
C
          RLM=ALM-TLM0
          SRLM= SIN(RLM)
          CRLM= COS(RLM)
          SAPH= SIN(APH)
          CAPH= COS(APH)
          CC=CRLM*CAPH
C
          ANUM=SRLM*CAPH
          DENOM=CTPH0*CC+STPH0*SAPH
C
          ALM= ATAN2(ANUM,DENOM)
          APH= ASIN(CTPH0*SAPH-STPH0*CC)
C
C-----IS ALM,APH INSIDE THE (ONE LINE REDUCED) MODEL DOMAIN?------------
C
      IF (ALM.LE.ALMWB .OR. ALM.GE.ALMEB .OR.
     &    APH.LE.APHSB .OR. APH.GE.APHNB)  THEN
             INSIDE=.FALSE.
             RETURN
      END IF
             INSIDE=.TRUE.
C
C
C---------X1,Y1 IS A COORDINATE SYSTEM WITH DLM,DPH AS LENGTH UNITS-----
      X1=(ALM-WB)*RDLM
      Y1=(APH-SB)*RDPH
C---------X2,Y2 ROTATED FOR +45 DEG. & TRANSLATED FOR IM1---------------
      X2=D50*( X1+Y1)
      Y2=D50*(-X1+Y1)+RIM1
C---------I2,J2 ARE COORDINATES OF SOUTHERN POINTS OF GRID BOXES--------
      I2=  INT(X2)
      J2=  INT(Y2)
C---------REMAINING PARAMETERS NEEDED TO KNOW THE POSITION--------------
C         OF THE TRANSFORMED POINT WITHIN THE MODEL REGION
      P=X2-I2
      Q=Y2-J2
C-----------------INDEX K CORRESPONDS TO I2,J2--------------------------
      JR=J2-IM1
      I3=I2-JR
      J3=I2+JR
      K=J3*IM-J3/2+(I3+2)/2
C-----------------------------------------------------------------------
                             RETURN
                             END
C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
