SUBROUTINE SETUP (domain_check_h, IPROJ, PHIC, XLONC, &
                                        TRUELAT1, TRUELAT2, &
#ifdef BKG
                        DIS, XCNTR, YCNTR, XN, POLE, PSI1,  C2)
#else
                        MAXNES, NESTIX, NESTJX, DIS, NUMC, NESTI, NESTJ, &
                        IXC, JXC, XCNTR, YCNTR, XN, POLE, PSI1,  C2, &
                        XIM11, XJM11) 
#endif
!------------------------------------------------------------------------------!
!
! MAP PROJECTION FACTORS CALCULATION
!
! INPUT:
! -----
!  IPROJ:     PROJECTION TYPE
!  PHIC:      CENTRAL LATITUDE
!  XLONC:     CENTRAL LONGITUDE
!  TRUELAT1:  TRUE LATITUDE 1
!  TRUELAT2:  TRUE LATITUDE 2
!  MAXNES:    NUMBER OF DOMAINS < 10
!  NESTIX (): HORIZONTAL DOMAINS DIMENSIONS IN I DIRECTION
!  NESTJX (): HORIZONTAL DOMAINS DIMENSIONS IN J DIRECTION
!  DIS ():    HORIZONTAL DOMAINS GRID SPACES IN METERS
!  NUMC ():   MOTHER DOMAIN ID (1= COARSE DOMAIN)
!  NESTI ():  DOMAINS' I-INDECES OF LOW LEFT CORNER IN MOTHER DOMAIN
!  NESTJ ():  DOMAINS' J-INDECES OF LOW LEFT CORNER IN MOTHER DOMAIN
!
! OUTPUT:
! ------
!  IXC:         COARSE DOMAIN GRID I DIMENSION
!  JXC:         COARSE DOMAIN GRID J DIMENSION
!  XCNTR:       X-COORDINATE OF COARSE DOMAIN CENTER
!  YCNTR:       Y-COORDINATE OF COARSE DOMAIN CENTER
!  XN:          CONE PROJECTION
!  POLE:        POLE LATITUDE
!  PSI1:       
!  C2:
!  XIM11 (10):  DOMAINS' INDECES OF LOW LEFT CORNER OF DOMAIN IDD
!  XJM11 (10):  DOMAINS' INDECES OF LOW LEFT CORNER OF DOMAIN IDD
!
! 
! Y.-R. GUO,       September 2000
! F. VANDENBERGHE, March 2001
!------------------------------------------------------------------------------!

      IMPLICIT NONE

      INTEGER, INTENT (inout) :: iproj
      REAL,    INTENT (inout) :: phic,  xlonc, truelat1, truelat2
#ifdef BKG
      REAL,    INTENT (out)   :: pole, xn, psi1, c2, xcntr, ycntr
      REAL,    INTENT (inout), DIMENSION (3) :: dis
#else
      INTEGER, INTENT (in)    :: maxnes
      INTEGER, INTENT (in), DIMENSION (maxnes) :: nestix, nestjx, NUMC, &
                                              nesti, nestj
      REAL,    INTENT (inout), DIMENSION (maxnes) :: dis

      INTEGER, INTENT (out)   :: ixc, jxc

      REAL,    INTENT (out)   :: pole, xn
      REAL,    INTENT (out)   :: xim11 (maxnes), xjm11 (maxnes)
      REAL,    INTENT (out)   :: psi1, c2, xcntr, ycntr
#endif
!------------------------------------------------------------------------------!
#ifdef BKG
      REAL :: psx, cell, cell2, r, phictr, sign

      INTEGER, DIMENSION (3) :: nproj

      CHARACTER (LEN=80) :: project
#else
!     REAL :: conv, a
      REAL :: iexp, aexp
      REAL :: psx, cell, cell2, r, phictr, hdsize, sign
      REAL :: xjc, xdis, yjc, ydis

      REAL :: cntrj0, cntri0
      INTEGER :: ixcmax, jxcmax
      REAL, DIMENSION (maxnes) :: xsouth,xwest,xnorth,xeast

      INTEGER :: incr
      INTEGER :: m, nm, nmc, nd1, nd2, mismatch
      INTEGER :: ioffst, joffst, iimxn, jjmxn
      INTEGER :: iratio (maxnes), nratio (maxnes)

      INTEGER, DIMENSION (3) :: nproj

      CHARACTER (LEN=80) :: project
!------------------------------------------------------------------------------!
!     DATA conv  / 57.29578 /
!     DATA a     / 6370.    /
#endif
      include 'constants.inc'
      
      logical  domain_check_h
      DATA nproj / 1, 2, 3  /
!     DATA nproj / 'LAMCON', 'POLSTR', 'MERCAT' /
!------------------------------------------------------------------------------!
      if (iproj == 0 ) then
          PROJECT='CE'
          XCNTR = 0
          YCNTR = 0
!          DIS(1) = 2.0 * pi * a / (NESTJX (1)-1) 
          DIS(1) = 110.0
          if (domain_check_h) STOP 'domain_check_h'
          goto 1100
      endif

!     WRITE (0,'(A)') " SET UP MAP PARAMETERS"

#ifdef BKG
! DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ:                            
                                                                               
      XN = -1.0E08                                                           

      IF (PHIC.LT.0) THEN                                                        
         SIGN=-1.      ! SOUTH HEMESPHERE 

      ELSE                                                                      

         SIGN=1.       ! NORTH HEMESPHERE

      ENDIF                                                                     

      POLE = 90.                                                             

      IF (IPROJ .EQ. NPROJ(1)) THEN
          if (abs(TRUELAT1 - TRUELAT2) .gt. 0.1) then
            XN = ALOG10 (COS (TRUELAT1 / CONV)) - &
                 ALOG10 (COS (TRUELAT2 / CONV))                                         
 
            XN = XN/(ALOG10 (TAN ((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - & 
                     ALOG10 (TAN ((45.0 - SIGN*TRUELAT2/2.0) / CONV)))          
          else
            XN = sign*sin(truelat1 / conv)
          endif

          PSI1 = 90.-SIGN*TRUELAT1

          PROJECT='LC'

      ELSE IF (IPROJ .EQ. NPROJ(2)) THEN

          XN = 1.0


! PSI1 IS THE PSEUDO-LATITUDE                                                   

          PSI1 = 90.-SIGN*TRUELAT1

          PROJECT = 'ST'

      ELSE IF (IPROJ .EQ. NPROJ(3)) THEN

          XN = 0.                                                                

          PSI1 = 0.
          PROJECT = 'ME'

      END IF                                                                    

      PSI1 = PSI1 / CONV

      IF (PHIC .LT. 0.) THEN                                                      
          PSI1 = -PSI1
          POLE = -POLE

      ENDIF                                                                     

!--------CALCULATE R                                                            

      IF (IPROJ .NE. NPROJ(3)) THEN

          PSX = (POLE-PHIC)/CONV

          IF (IPROJ .EQ. NPROJ(1)) THEN

              CELL  = A*SIN (PSI1)/XN
              CELL2 = (TAN (PSX/2.)) / (TAN (PSI1/2.))                               
         ENDIF                                                                  

         IF (IPROJ .EQ. NPROJ(2)) THEN
            CELL  = A*SIN (PSX)/XN
            CELL2 = (1. + COS (PSI1))/(1. + COS (PSX))
         ENDIF                                                                  

         R = CELL*(CELL2)**XN                                                   

         XCNTR = 0.0                                                            
         YCNTR = -R                                                             

      ENDIF                                                                     

!-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1            

      IF (IPROJ .EQ. NPROJ(3)) THEN

         C2     = A*COS (PSI1)
         XCNTR  = 0.0                                                           
         PHICTR = PHIC/CONV                                                     
         CELL   = COS (PHICTR)/(1.0+SIN (PHICTR)) 
         YCNTR  = - C2*ALOG (CELL)                                              

      ENDIF                                                                     
1100  continue

#else
                                                                              
! REMOVE DOMAIN EXTENSION IN MM5 INPUT/OUTPUT FILES

      IEXP = 0
      AEXP =  -1.001*DIS(1)

! REMOVE INITIALIZATION OF TRUE LAT 

!     TRUELAT1=91.0                                                             
!     TRUELAT2=91.0                                                             

! DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ:                            
                                                                               
      XN = -1.0E08                                                           

      IF (PHIC.LT.0) THEN                                                        
         SIGN=-1.      ! SOUTH HEMESPHERE 

      ELSE                                                                      

         SIGN=1.       ! NORTH HEMESPHERE

      ENDIF                                                                     

      POLE = 90.                                                             

      IF (IPROJ .EQ. NPROJ(1)) THEN

          IF (ABS (TRUELAT1) .GT. 90.) THEN
              TRUELAT1 = 60.
              TRUELAT2 = 30.
              TRUELAT1 = SIGN*TRUELAT1
              TRUELAT2 = SIGN*TRUELAT2
          ENDIF

          IF (ABS(TRUELAT1 - TRUELAT2) .GT. 0.1) then
            XN = ALOG10 (COS (TRUELAT1 / CONV)) - &
                 ALOG10 (COS (TRUELAT2 / CONV))

            XN = XN/(ALOG10 (TAN ((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - &
                     ALOG10 (TAN ((45.0 - SIGN*TRUELAT2/2.0) / CONV)))
          ELSE
            XN = SIGN*SIN(TRUELAT1 / CONV)
          ENDIF

          PSI1 = 90.-SIGN*TRUELAT1

          PROJECT='LC'

      ELSE IF (IPROJ .EQ. NPROJ(2)) THEN

          XN = 1.0

          IF (ABS(TRUELAT1) .GT. 90.) THEN

              TRUELAT1 = 60.
              TRUELAT2 = 0.

              TRUELAT1 = SIGN*TRUELAT1
              TRUELAT2 = SIGN*TRUELAT2

          ENDIF


! PSI1 IS THE PSEUDO-LATITUDE                                                   

          PSI1 = 90.-SIGN*TRUELAT1

          PROJECT = 'ST'

      ELSE IF (IPROJ .EQ. NPROJ(3)) THEN

          XN = 0.                                                                

          IF (ABS (TRUELAT1) .GT. 90.) THEN

              TRUELAT1 = 0.
              TRUELAT2 = 0.                                                               
          ENDIF

          IF (TRUELAT1 .NE. 0.) THEN                                                   
              WRITE (0,'(/,A)') &
             "MERCATOR PROJECTION IS ONLY SUPPORTED AT 0 DEGREE TRUE LATITUDE."
              WRITE (0,'(A,/)') &
             "TRUELAT1 IS RESET TO 0"

               TRUELAT1 = 0.                                                               
          ENDIF

          PSI1 = 0.
          PROJECT = 'ME'

      END IF                                                                    

      PSI1 = PSI1 / CONV

      IF (PHIC .LT. 0.) THEN                                                      
          PSI1 = -PSI1
          POLE = -POLE

      ENDIF                                                                     

!--------CALCULATE R                                                            

      IF (IPROJ .NE. NPROJ(3)) THEN

          PSX = (POLE-PHIC)/CONV

          IF (IPROJ .EQ. NPROJ(1)) THEN

              CELL  = A*SIN (PSI1)/XN
              CELL2 = (TAN (PSX/2.)) / (TAN (PSI1/2.))                               
         ENDIF                                                                  

         IF (IPROJ .EQ. NPROJ(2)) THEN
            CELL  = A*SIN (PSX)/XN
            CELL2 = (1. + COS (PSI1))/(1. + COS (PSX))
         ENDIF                                                                  

         R = CELL*(CELL2)**XN                                                   

         XCNTR = 0.0                                                            
         YCNTR = -R                                                             

      ENDIF                                                                     

!-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1            

      IF (IPROJ .EQ. NPROJ(3)) THEN

         C2     = A*COS (PSI1)
         XCNTR  = 0.0                                                           
         PHICTR = PHIC/CONV                                                     
         CELL   = COS (PHICTR)/(1.0+SIN (PHICTR)) 
         YCNTR  = - C2*ALOG (CELL)                                              

      ENDIF                                                                     

1100  continue

      WRITE (0,'(2(A,F8.1),A,2X,A,f10.3)') &
      " COARSE GRID CENTER IS AT X = ",XCNTR," KM AND Y = ",YCNTR," KM.", &
      " DIS(1)=", DIS(1) 


!   CHECK THE COMPATIBILITY OF NEST DOMAINS WITH THE COARSE DOMAINS             
!     AND CALCULATE THE IRATIOS, INORTHS, ISOUTHS, IWESTS AND IEASTS            

!     A) EXTENDING THE COARSE DOMAIN IF IEXP = 1                                

        IXC = NESTIX (1)
        JXC = NESTJX (1)

        IOFFST = 0                                                              
        JOFFST = 0                                                              

      IF (IEXP .EQ. 1) THEN

          INCR = INT (AEXP/DIS (1) + 1.001)

          IXC = NESTIX (1) + INCR*2
          JXC = NESTJX (1) + INCR*2

          IOFFST = INCR
          JOFFST = INCR

          WRITE (0,'(A,I3)') " GRID IS EXPANDED BY ",INCR, &
                             " GRID POINTS ON EACH EDGE."
!          WRITE (0,'(2(A,I3))') &
!         "EXTENDED IXC IS ",IXC," EXTENDED JXC IS ",JXC               

      ENDIF                                                                     

!-----CENTER OF GRID IN THE COARSE DOMAIN                                       

      CNTRJ0 = FLOAT (JXC+1)/2.
      CNTRI0 = FLOAT (IXC+1)/2.

!      WRITE (0,'(2(A,I5))') &
!     "COARSE DOMAIN SIZE IX = ",NESTIX(1)," JX = ", NESTJX(1)                                            
!  MIX, MJX ARE USED IN SUB. TFUDGE:                                            

      IXCMAX = IXC
      JXCMAX = JXC 

      DO M = 1, MAXNES                                                       
         IXCMAX = MAX0 (NESTIX(M),IXCMAX)
         JXCMAX = MAX0 (NESTJX(M),JXCMAX)
      ENDDO

!      PRINT 24, IXCMAX,JXCMAX
!24    FORMAT('   ++> THE MAXIMUM DIMENSION = (',2I5,') <++')

!  CHECK IF POLE IS INSIDE THE DOMAIN OR NOT FOR LAMBERT PROJECTION:            

      HDSIZE = (IXC-1)*DIS(1)/2.                                               

!     IF (HDSIZE.GT.ABS(YCNTR) .AND. IPROJ.EQ.NPROJ(1))THEN                     
!        PRINT *,'-------------------------------------------------'            
!        PRINT *,'HALF DOMAIN SIZE IN Y-DIRECTION = ',HDSIZE                    
!        PRINT *,'    DISTANCE FROM CNTER TO POLE = ',ABS(YCNTR)                
!        PRINT *,'NOT MAKE SENSE WITH THE POLE IS INSIDE THE DOMAIN '           
!        PRINT *,'    FOR LAMBERT CONFORMAL PROJECTION!'                        
!        PRINT *,'=== PLEASE RE-SPECIFY THE CENTER OR DOMAIN SIZE. ==='         
!        STOP                                                                   
!     ENDIF                                                                     

!     B) CALCULATING THE IRATIOS, INORTHS AND IEASTS:                           

      IRATIO (1) = 1
      NRATIO (1) = 1
      XSOUTH (1) = 1.
      XWEST  (1) = 1.
      XNORTH (1) = FLOAT (IXC)
      XEAST  (1) = FLOAT (JXC)

      XJC = (XEAST(1) + 1.0)/2.                                                 

!      PRINT 27,XSOUTH(1),XWEST(1),XNORTH(1),XEAST(1),DIS(1),& 
!               IRATIO(1),NRATIO(1)                                            
! 27   FORMAT(1X,'XSOUTH(1)= ',F6.1,2X,'XWEST(1)= ',F6.1,2X, &
!      'XNORTH(1)= ',F6.1,2X,'XEAST(1)= ',F6.1,2X,'DIS(1)= ',F6.1,&
!      2X,'IRATIO(1)= ',I3,2X,'NRATIO(1)= ',I3)                                  

      MISMATCH = 0                                                              
                                                                                
      DO 30 NM = 2, MAXNES                                                      

!  DOMAINS' CONSISTENCY CHECK:                                                  

      NMC = NUMC (NM)

      IF (AMOD ((DIS (NMC)+0.09),DIS (NM)) .GT. 0.1) THEN

         MISMATCH = MISMATCH + 1                                                

!        PRINT 29, NM,NMC                                                       
!        PRINT 31,NM,DIS(NM),NMC,DIS(NMC)                                       
!29      FORMAT(2X,'DOMAIN ',I2,' HAS INCORRECT GRID SIZE ', &
!      '  IT HAS TO BE THE MULTIPLE OF DOMAIN ',I2)                             
!31      FORMAT(2X,'DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM', &
!               '  DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM')                        

        GO TO 30                                                                

      ENDIF                                                                     

      IRATIO (NM) = NINT (DIS (NMC)/ DIS (NM))
      NRATIO (NM) = NINT (DIS (1)  / DIS (NM))                                         

!   MAKE SURE THE 4 CORNER POINTS OF NEST DOMAINS ARE ON THE                    
!   PREVIOUS DOMAIN GRID-POINTS                                                 

      IF (MOD((NESTIX(NM)-1),IRATIO(NM)).NE.0) THEN                             

        MISMATCH = MISMATCH + 1                                                 

        IIMXN = (INT(FLOAT(NESTIX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1          

!        PRINT 32,NM,NESTIX(NM),IRATIO(NM),IIMXN                                
!32      FORMAT(2X,'NESTIX(',I2,')=',I4,' AND IRATIO=',I2,&
!      ' DOES NOT MATCH, YOU MAY SET NESTIX TO ',I4)

      ENDIF                                                                     

      IF (MOD ((NESTJX (NM)-1),IRATIO (NM)).NE.0) THEN

         MISMATCH = MISMATCH + 1

        JJMXN = (INT(FLOAT(NESTJX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1          

!        PRINT 33,NM,NESTJX(NM),IRATIO(NM),JJMXN                               
!33      FORMAT(2X,'NESTJX(',I2,')=',I4,' AND IRATIO=',I2,&
!      ' DOES NOT MATCH, YOU MAY SET NESTJX TO ',I4)                                         
      ENDIF                                                                    
!                                                                               
!-----REDEFINE LOCATION OF LOWER LEFT CORNER OF FINE MESH (IN TERMS          
!     OF EXTENDED COARSE MESH - DOMAIN 1 INDICES) IF USING EXTENDED G        
                                                                                
!      PRINT 34,NM,NESTIX(NM),NESTJX(NM),DIS(NM),NESTI(NM),NESTJ(NM),&
!      NUMC(NM),IRATIO(NM),NRATIO(NM),IEXP                                   
!34    FORMAT(/1X,'DOMAIN ',I2,2X,'IX=',I3,2X,'JX=',I3,2X, &
!      'DS= ',F6.1,2X,'ICNS=',I4,2X,'JCNS=',I4,2X,&
!      'NUMC= ',I2,2X,'IRATIO= ',I2,2X,'NRATIO= ',I2,2X,&
!      'IEXP= ',I2)                                                   

      XDIS = 0.0                                                                
      YDIS = 0.0                                                                
      ND1 = NM                                                                  
      ND2 = NMC                                                                 

40    CONTINUE                                                                  

      XDIS = (NESTI(ND1)-1)*DIS(ND2) + XDIS
      YDIS = (NESTJ(ND1)-1)*DIS(ND2) + YDIS

      IF (ND2 .GT. 1) THEN                                                      
        ND1 = ND2                                                               
        ND2 = NUMC(ND2)                                                        
        GO TO 40                                                                
      ENDIF                                                                     

      XSOUTH(NM) = XDIS/DIS(1) + FLOAT(IOFFST) + 1                             
      XWEST(NM)  = YDIS/DIS(1) + FLOAT(JOFFST) + 1                             
      XNORTH(NM) = XSOUTH(NM) + FLOAT(NESTIX(NM)-1)*DIS(NM)/DIS(1)             
      XEAST(NM)  = XWEST(NM) + FLOAT(NESTJX(NM)-1)*DIS(NM)/DIS(1)             
                                                                                
!      PRINT 35                                                                 
!      PRINT 36,XSOUTH(NM),XWEST(NM),XNORTH(NM),XEAST(NM)                       
!35    FORMAT(2X,'COARSE MESH INDICES FOR THE 4 CORNER POINTS ARE')             
!36    FORMAT(2X,'SOUTH=',F6.1,3X,'WEST=',F6.1,3X,'NORTH =',F6.1,3X, & 
!            'EAST = ',F6.1)                                                   

30    CONTINUE                                                                 

      IF (MISMATCH.GT.0) THEN                                                  

!       PRINT *, &
!      'TERRAIN STOP IN SUBROUTINE SETUP DUE TO INCORRECT NEST DOMAIN SET UP' 
!       STOP 1111

       ENDIF                                                                    

!      Copy output

       
       XIM11 = XSOUTH
       XJM11 = XWEST

#endif
       END SUBROUTINE SETUP