!
!***********************************************************************
!                                                                      *
   SUBROUTINE SWREAD (COMPUT)                                          
!                                                                      *
!***********************************************************************
!
!     Reading and processing of the user commands describing the model
!
!***********************************************************************

   USE TIMECOMM                                                        
   USE OCPCOMM1                                                        
   USE OCPCOMM2                                                        
   USE OCPCOMM4                                                        
   USE SWCOMM1                                                         
   USE SWCOMM2                                                         
   USE SWCOMM3                                                         
   USE SWCOMM4                                                         
   USE OUTP_DATA                                                       
   USE M_SNL4                                                          !
   USE M_GENARR                                                        
   USE M_OBSTA                                                         
!   USE M_PARALL
   USE MOD_UTILS
!   USE MOD_USGRID
#  if defined (EXPLICIT)
   USE MOD_ACTION_EX
#  else      
   USE MOD_ACTION_IM
#  endif   
#  if defined (MULTIPROCESSOR)   
   USE MOD_PAR
#  endif  
   USE ALL_VARS
   USE VARS_WAVE
   USE MOD_OBCS
   USE MOD_NCLL
   USE MOD_NCTOOLS
   USE MOD_FORCE

   IMPLICIT NONE
!
!     SWAN (Simulating WAves Nearshore); a third generation wave model
!     Copyright (C) 2008  Delft University of Technology
!     FVCOM-SWAVE; a third generation wave model
!     Copyright (C) 2008-2009  University of Massachusetts Dartmouth
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation; either version 2 of
!     the License, or (at your option) any later version.
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     A copy of the GNU General Public License is available at
!     http://www.gnu.org/copyleft/gpl.html#SEC3
!     or by writing to the Free Software Foundation, Inc.,
!     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
!
!
!  0. Authors
!     Jianhua Qi
!
!  1. Updates
!
!  2. Purpose
!
!     Reading and processing of the user commands describing the model
!
!  3. Method (updated 30.72)
!
!  4. Argument variables (updated 30.72)
!
!  5. Parameter variables
!
!  6. Local variables
!

   INTEGER, PARAMETER :: MXINCL = 10                                   

   INTEGER, SAVE     :: INCNUM(1:MXINCL) = 0                           
   INTEGER, SAVE     :: INCLEV = 1                                     
!
   INTEGER           :: IERR = 0                                       
   INTEGER           :: ICNL4, ILAMBDA                                 
   INTEGER           :: ITMP1, ITMP2, ITMP3, ITMP4                     
!
   LOGICAL           :: MORE                                           
   LOGICAL, SAVE     :: LOBST = .FALSE.                                
!
   INTEGER        :: LREF, LREFDIFF, LRFRD                             
   REAL           :: KDIF                                              
   REAL           :: POWD                                              
   REAL           :: POWS, DUM                                         
   REAL           :: FD1, FD2, FD3, FD4                                
   REAL, ALLOCATABLE :: RLAMBDA(:)                                     

   TYPE(OBSTDAT), POINTER :: OBSTMP                                    
   TYPE(OBSTDAT), SAVE, POINTER :: COBST                               

   TYPE(OPSDAT), POINTER :: OPSTMP                                     

   TYPE XYPT                                                           
     REAL                :: X, Y
     TYPE(XYPT), POINTER :: NEXTXY
   END TYPE XYPT

   TYPE(XYPT), TARGET  :: FRST                                         
   TYPE(XYPT), POINTER :: CURR, TMP                                    

   LOGICAL :: KEYWIS                                                   
   LOGICAL :: EQREAL                                                   

   LOGICAL :: FOUND                                                    
   LOGICAL, SAVE :: RUNMADE = .FALSE.                                  

   LOGICAL, SAVE :: LOGCOM(1:6) = .FALSE.                              
   INTEGER   TIMARR(6)                                                 
   CHARACTER PSNAME *8, PNAME *8, COMPUT *(*), PTYPE *1, DTTIWR *18    
   INTEGER, SAVE :: IENT = 0       ! number of entries to this subr
   INTEGER, SAVE :: LWINDR = 0     ! if non-zero, there is wind
   INTEGER, SAVE :: LWINDM = 3     ! type of wind growth formulation
   INTEGER, ALLOCATABLE :: IARR(:) 
   INTEGER :: GEN,IOS,LMAX                                    

   CHARACTER(LEN=7) INPFIL
   CHARACTER(LEN=80) CHTMP
   CHARACTER(LEN=80) PROJJ,SET,MODE,COORDINATES,BOUNDARY,CONSTANT,FRIC_FORM
   REAL FNTMP
   INTEGER INTMP
   LOGICAL LTMP,STAT,CART,CHECK
   INTEGER  INTVEC(150),ISCAN   !,KTEMP
   
   CHARACTER(LEN=80) INP_BOT_NAME,INP_CUR_NAME,INP_WI_NAME,      &
                     INP_FR_NAME,INP_WLEV_NAME,PROP_CHOICE
		     
   CHARACTER(LEN=20) CHARWAVEPERIOD
   CHARACTER(LEN=20) GROWTH
   CHARACTER(LEN=20) CNL4,OBST_TYPE,RDIFF,OFF_TYPE,NUM_CHOICE,STATL, &
                     SIGIM
   CHARACTER(LEN=20) WCAP		     
   CHARACTER(LEN=6)  UNIT
   LOGICAL QUAD,AGROW,BRE,MDIAL,FRICL,OBST,TRI,REFL,RFD,SETUP,    &
           DIFFR,OFF,PROP,FLUXLIM,DIR,COM_STAT	
   LOGICAL INP_BOT,INP_CUR,INP_CUR_STAT,INP_CUR_SERI,INP_WI,INP_WI_STAT, &
           INP_WI_SERI,INP_FR,INP_FR_STAT,INP_FR_SERI,INP_WLEV,          &
	   INP_WLEV_STAT,INP_WLEV_SERI
   LOGICAL LAM,LIM,LINE,NUM	   
   	   
	   
!   REAL, ALLOCATABLE :: DEPTH_GL(:),UXB_GL(:,:),UYB_GL(:,:),             &
!           WXI_GL(:,:),WYI_GL(:,:),FRIC_GL(:,:),WLEVL_GL(:,:)
   REAL, ALLOCATABLE :: UXB_GL(:,:),UYB_GL(:,:),FRIC_GL(:,:),WLEVL_GL(:,:)
   REAL :: IFLTMP,NTIME_TMP,TEMPXY	   
   
   REAL :: ALTMP,FRLOW,FRHIG,GAMMA,TMPDIR,TRCF,HGT,OGAM,OBET,REF0,DIFF
   INTEGER :: MSS,IVAL
   INTEGER :: I,IGR2,IGRD,IGR1,J,JJ,NUMCOR,I1,NCNT
   INTEGER :: ITRAS	   	      
   INTEGER, ALLOCATABLE :: TEMP1(:),TEMP2(:)
   REAL    :: DEGCNV
!
!  ***** read command *****
!
!
!  ------------------------------------------------------------------
!
!  PROJECT      reading of project title and description
!
! ===============================================================
!
! PROJect  'NAME'  'NR'
!
!          'title1'
!
!          'title2'
!
!          'title3'
!
! ===============================================================
!

   CALL INCSTR('PROJECT',PROJJ,'REQ',' ')
   
   IF(PROJJ /= 'default')THEN
     CALL INCSTR('NAME',PROJID,'UNC',' ')
     CALL INCSTR('NR',PROJNR,'REQ',' ')
     CALL INCSTR('title1',PROJT1,'UNC',' ')   
     CALL INCSTR('title2',PROJT2,'UNC',' ')    
     CALL INCSTR('title3',PROJT3,'UNC',' ')
   END IF  
    
!
!  ------------------------------------------------------------------
!
!  SET       setting physical parameters and error counters
!
! =============================================================
!
!   SET  [level]  [nor]  [depmin]  [maxmes]        &
!        [maxerr]  [grav]  [rho]  [inrhog]         &
!        [hsrerr]  CARTesian/NAUTical  [pwtail]    &
!        [froudmax]  [printf]  [prtest]
!
! =============================================================
!
   CALL INCSTR('SET',SET,'REQ',' ')
   
   IF(SET /= 'default')THEN
     CALL INREAL('LEVEL', WLEV,  'UNC',0.0)
     CALL INREAL('NOR',   DNORTH,'UNC',0.0)
     CALL INREAL('DEPMIN',DEPMIN,'UNC',0.0)
     CALL ININTG('MAXMES',MAXMES,'UNC',0)
     CALL ININTG('MAXERR',MAXERR,'UNC',0)
     CALL INREAL('GRAV',  GRAV_W,  'UNC',0.0)
     CALL INREAL('RHO',   RHO_W,   'UNC',0.0)
     CALL ININTG('INRHOG',INRHOG,'UNC',0)
     IF(INRHOG == 0)THEN                                            
       OVUNIT(7) = 'm2/s'
       OVUNIT(9) = 'm2/s'
       OVUNIT(19) = 'm3/s'
       OVUNIT(21) = 'm2s'
       OVUNIT(22) = 'm2'
       OVUNIT(29) = 'm2'
     ELSE
       OVUNIT(7) = 'W/m2'
       OVUNIT(9) = 'W/m2'
       OVUNIT(19) = 'W/m'
       OVUNIT(21) = 'Js/m2'
       OVUNIT(22) = 'J/m2'
       OVUNIT(29) = 'J/m2'
     ENDIF                                                            
      
     CALL INREAL('HSRERR',HSRERR,'UNC',0.0)

     CALL INLOGC('NAUTICAL',BNAUT,'UNC','T')
     IF(BNAUT) OUTPAR(4) = 0.   
     IF(ITEST >= 20)THEN                                        
       WRITE(PRTEST,*) ' Set NAUT command: BNAUT=', BNAUT            
     END IF                                                            

     CALL INREAL('PWTAIL',PWTAIL(1),'UNC',0.0)
     IF(CHGVAL)THEN
       IF(PWTAIL(1) <= 1.) CALL MSGERR (3, 'Incorrect PWTAIL')
       PWTAIL(3) = PWTAIL(1) + 1.
     END IF
     
     CALL INREAL('FROUDMAX',PNUMS(18),'UNC',0.0)
     CALL ININTG('PRINTF',  PRINTF,   'UNC',0)
     CALL ININTG('PRTEST',  PRTEST,   'UNC',0)
   END IF 
   
!
!  ------------------------------------------------------------------
!
!  MODE  : Set STATionary, DYNamic (NONSTAtionary) or 1D SWAN model    
!
! ========================================
!
!        | -> STAtionary    |     | -> TWODimensional |                   
! MODE  <                    >   <                     >  (NOUPDATe)      
!        |    NONSTationary |     |    ONEDimensional |                   
!
! ========================================
!
   CALL INCSTR('MODE',MODE,'REQ',' ')

   IF(MODE /= 'default')THEN
     CALL INLOGC('STATIONARY',STAT,'UNC','T')

     IF(.NOT. STAT)THEN
       IF(NSTATM == 0) CALL MSGERR (2, 'Mode Nonst incorrect here')   
       NSTATM = 1                                                      
       NSTATC = 1                                                      
!
!      switch on flag for computation of default initial condition     
       ICOND = 1                                                       
!      no relaxation                                                   
       PNUMS(30) = 0.                                                  
     ELSE  
       IF(NSTATM == 1) CALL MSGERR (2, 'Mode STAT incorrect here')    
       NSTATM = 0                                                      
     ENDIF                                                             

     CALL INLOGC('ONED',ONED,'UNC','T')
   END IF

   CALL INLOGC('ACUPDAT',ACUPDA,'UNC','F')

!
!  ------------------------------------------------------------------
!
!  COORD : spherical or cartesian coordinates
!
! =============================================================
!
!   COORDinates  /  -> CARTesian               \   REPeating
!                \ SPHErical [rearth]   UM/QC  /
!
! =============================================================
!
   CALL INCSTR('COORDINATES',COORDINATES,'REQ',' ')
   
   IF(COORDINATES /= 'default')THEN
     CALL ININTG('KSPHER',KSPHER,'UNC',0)

     IF(KSPHER ==1)THEN
#    if defined (SPHERICAL)
       CALL INREAL ('REARTH', REARTH2, 'UNC', 0.)  !Jianzhong Ge                     
       LENDEG = REARTH2 * PI_W / 180.                                   
!      change properties of output quantities Xp and Yp
       OVUNIT(1) = 'degr'
       OVLLIM(1) = -200.
       OVULIM(1) =  400.
       OVLEXP(1) = -180.
       OVHEXP(1) =  360.
       OVEXCV(1) = -999.
       OVUNIT(2) = 'degr'
       OVLLIM(2) = -100.
       OVULIM(2) =  100.
       OVLEXP(2) = -90.
       OVHEXP(2) =  90.
       OVEXCV(2) = -999.

       CALL ININTG('PROJ_METHOD',PROJ_METHOD,'UNC',0)
#      endif
     ELSE IF(KSPHER /= 0)THEN
       WRITE(PRINTF,*)'PARAMETER KSPHER NOT EQUAL TO 0 OR 1: ',KSPHER
       CALL PSTOP
     END IF
   END IF
   
   CALL ININTG('KREPTX',KREPTX,'STA',0)
       
!
!     ------------------------------------------------------------------
!
!     CGRID     definition of computational grid
!
! ===========================================================================
!
!   CGRID  / REGular [xpc] [ypc] [alpc] [xlenc] [ylenc] [mxc] [myc] \
!          \ CURVilinear [mxc] [myc]  [excval]                      /    &
!          / CIRcle               \
!          \ SECtor [dir1] [dir2] /  [mdc]  [flow]  [fhig]  [msc]
!
! ===========================================================================
!
!   CALL INCSTR('GRIDFILE',GRIDFILE,'REQ',' ')
!   INQUIRE(FILE=TRIM(GRIDFILE),EXIST=CHECK) 
!   IF(CHECK)THEN
!     OPEN(1, FILE=TRIM(GRIDFILE),STATUS='OLD')
!   
!   ELSE
!     WRITE(PRINTF,*) TRIM(GRIDFILE), ' DOES NOT EXIST.'
!     PRINT*, TRIM(GRIDFILE), ' DOES NOT EXIST.'
!     CALL PSTOP
!   END IF
   
   LOGCOM(2) = .TRUE.   
   
   OPTG = 1

   CALL INREAL('ALPC',ALPC,'UNC',0.0)
!
!  ALPC is made to be between -PI and PI
!
   ALTMP = ALPC / 360.
   ALPC = PI2_W * (ALTMP - NINT(ALTMP))
   CVLEFT = .TRUE.                                                 

   CALL INLOGC('FULCIR',FULCIR,'UNC','T')
   IF(.NOT. FULCIR)THEN
     CALL INREAL('DIR1',SPDIR1,'REQ',0.0)
     CALL INREAL('DIR2',SPDIR2,'REQ',0.0)
   END IF  

   CALL ININTG('MDC',MDC,'REQ',0)
   CALL INREAL('FLOW',FRLOW,'REQ',0.0)
   CALL INREAL('FHIGH',FRHIG,'REQ',0.0)
   CALL ININTG('MSC',MSS,'REQ',0)
   IVAL = MAX(ABS(MSS*NINT(FRLOW)),ABS(MSS*NINT(FRHIG)),              &
              ABS(NINT(FRLOW)*NINT(FRHIG)))                          
   IF(IVAL >= 999**2)THEN                                          
     CALL MSGERR(3,'At least, FLOW and FHIGH or MSC must be given!')  
   END IF                                                            
   IVAL = MAX(ABS(NINT(FRLOW)),ABS(NINT(FRHIG)),ABS(MSS))            
   IF(IVAL /= 999)THEN                                             
     GAMMA = EXP(ALOG(FRHIG/FRLOW)/REAL(MSS)) - 1.                  
     WRITE (PRINTF,'(A,(F7.4))')                                      &
                  ' Resolution in sigma-space: df/f = ',GAMMA  
   ELSE IF(MSS == -999)THEN                                        
     MSS = NINT(ALOG(FRHIG/FRLOW)/ALOG(1.1))                        
     WRITE (PRINTF,'(A,I3)')                                          &
                  ' Number of meshes in sigma-space: MSC-1 = ',MSS  
   ELSE IF(FRLOW == -999.)THEN                                     
     FRLOW = FRHIG/EXP(ALOG(1.1)*REAL(MSS))                         
     WRITE (PRINTF,'(A,(F7.4))')                                      &
                  ' Lowest discrete frequency (in Hz): FLOW = ',FRLOW  
   ELSE IF(FRHIG == -999.)THEN                                     
     FRHIG = FRLOW*EXP(ALOG(1.1)*REAL(MSS))                         
     WRITE (PRINTF,'(A,(F7.4))')                                      &
                  ' Highest discrete frequency (in Hz): FHIGH = ',FRHIG  
   END IF                                                            
   SLOW = 2.*PI_W*FRLOW
   SHIG = 2.*PI_W*FRHIG
   MSC  = MSS+1
!
!    ***** MCGRD is the total number of nodes   *****
!    ***** MDC is the number of steps in Theta direction
!              as parts of a circle
!    ***** MSC is the number of frequencies in sigma direction (logaritmic)
     IF(FULCIR)THEN
       DDIR  = PI2_W / MDC                                               
!
!      modification of SPDIR1 first installed with version 30.72, then
!      reversed and reinstalled with 40.13
!      purpose: prevent problem with grids under 45 degrees            
       SPDIR1 = ALPC + 0.5 * DDIR                                      
     ELSE
       IF(BNAUT)THEN                                                 
!        swap values of SPDIR1 and SPDIR2, and transform               
         TMPDIR = SPDIR1                                               
         SPDIR1 = 180. + DNORTH - SPDIR2                               
         SPDIR2 = 180. + DNORTH - TMPDIR                               
       ENDIF                                                           
       SPDIR1 = SPDIR1 * PI_W / 180.                                     
       SPDIR2 = SPDIR2 * PI_W / 180.                                     
       IF(SPDIR2 < SPDIR1) SPDIR2 = SPDIR2 + PI2_W                     
       DDIR = (SPDIR2-SPDIR1) / REAL(MDC)                              
       MDC = MDC + 1                                                   
     ENDIF
!
     IF(.NOT.ALLOCATED(SPCSIG)) ALLOCATE(SPCSIG(MSC))                  
     IF(.NOT.ALLOCATED(SPCDIR)) ALLOCATE(SPCDIR(MDC,6))                
     CALL SSFILL(SPCSIG,SPCDIR)                                        
!
     IF(ITEST >= 20)THEN
       IF(OPTG == 1)WRITE (PRINTF,6048)                              
       IF(OPTG == 3)WRITE (PRINTF,6049)                              
       WRITE (PRINTF,6045) SLOW, SHIG, FRINTF
       WRITE (PRINTF,6046) MCGRD,NGL,MDC,MSC
       WRITE (PRINTF,6047) DDIR
!       WRITE (PRINTF,6047) DX,DY,DDIR
6048   FORMAT ('GRID: REGULAR RECTANGULAR')                            
6049   FORMAT ('GRID: CURVILINEAR')                                    
6045   FORMAT (' S-low: ', F6.3,' S-hig: ', F6.3, ' frintf: ', F6.3)
6046   FORMAT (' MCGRD: ',I6,' MDC: ',I6,' MSC: ',I6)
6047   FORMAT (' DDIR: ', F6.3)
     ENDIF
!
     IF(.NOT.ALLOCATED(XCGRID)) ALLOCATE(XCGRID(MT))              
     IF(.NOT.ALLOCATED(YCGRID)) ALLOCATE(YCGRID(MT))              
!
     IF(OPTG == 1)THEN                                             
!      *** The coordinates of computational points in ***
!      *** regular grid are computed                  ***
       COSPC = COS(ALPC)                                               
       SINPC = SIN(ALPC)                                               
       DO I = 1,MT
         XCGRID(I)=VX(I)                                            
         YCGRID(I)=VY(I)                                            
       END DO   
     ENDIF
!
!    *** the computational grid is included in output data  ***        
     IF(MCGRD > 0)THEN                            
       ALLOCATE(OPSTMP)                                             
       OPSTMP%PSNAME = 'COMPGRID'                                   
       IF(OPTG == 1)THEN                                          
         OPSTMP%PSTYPE = 'F'                                        
         OPSTMP%OPR(1) = VXMIN       !XPC                                        
         OPSTMP%OPR(2) = VYMIN       !YPC                                        
         OPSTMP%OPR(3) = 1.0         !XCLEN                                      
         OPSTMP%OPR(4) = 1.0         !YCLEN                                      
         OPSTMP%OPR(5) = ALPC                                       
       ELSE                                                         
         OPSTMP%PSTYPE = 'H'                                        
         OPSTMP%OPR(1) = FLOAT(MXC-1)                                
         OPSTMP%OPR(2) = FLOAT(MYC-1)                                
         OPSTMP%OPR(3) = 0.                                          
         OPSTMP%OPR(4) = 0.                                          
         OPSTMP%OPR(5) = ALPC                                        
       ENDIF                                                          
       OPSTMP%OPI(1) = MCGRD       !MXC                                            
       ALLOCATE(OPSTMP%XP(0))                                         
       ALLOCATE(OPSTMP%YP(0))                                         
       NULLIFY(OPSTMP%NEXTOPS)                                        
       IF(.NOT.LOPS)THEN                                          
         FOPS = OPSTMP                                               
         COPS => FOPS                                                
         LOPS = .TRUE.                                               
       ELSE                                                           
         COPS%NEXTOPS => OPSTMP                                      
         COPS => OPSTMP                                              
       END IF                                                         
     ENDIF
     
!     ------------------------------------------------------------------
!
!     INPUT    definition of input grids
!
! ============================================================================
!
!   INPgrid                                                                &
!      BOTtom / WLEVel / CURrent / VX / VY / FRiction / WInd / WX / WY     &
!      / REG [xpinp] [ypinp] [alpinp]  [mxinp] [myinp]  [dxinp] [dyinp] \
!      \ CURVilinear [stagrx] [stagry] [mxinp] [myinp]                  /  &
!      (NONSTATionary [tbeginp] [deltinp] SEC/MIN/HR/DAY [tendinp])
!
! ============================================================================
   IGR2 = 0
   INP_BOT = .TRUE.
   IF(INP_BOT)THEN
     IGRD = 1
     IGR1 = 1
     LOGCOM(3) =.TRUE.                                                 
     PSNAME = 'BOTTGRID'
    
     ALLOCATE(DEPTH(MT))      ;  DEPTH    = 0.0

     DEPTH = H
   END IF

# if defined (WAVE_CURRENT_INTERACTION) && !defined (OFFLINE)
    ICUR = 1
# endif 
   CALL INLOGC('INP_CUR',INP_CUR,'REQ','T')
   IF(INP_CUR)THEN
     IGRD = 1
     IGR1  = 2
     IGR2 = 3
     ICUR   = 1
     PSNAME = 'VXGRID  '

     CALL INCSTR('INP_CUR_NAME',INP_CUR_NAME,'REQ',' ')
     INQUIRE(FILE=TRIM(INP_CUR_NAME),EXIST=CHECK)
     IF(CHECK)THEN
       OPEN(1,FILE=TRIM(INP_CUR_NAME),STATUS='OLD',FORM='UNFORMATTED')
     ELSE
       WRITE(PRINTF,*) TRIM(INP_CUR_NAME), ' DOES NOT EXIT.'
       CALL PSTOP
     END IF
    
     CALL INLOGC('INP_CUR_STAT',INP_CUR_STAT,'REQ','T')
     IF(.NOT. INP_CUR_STAT)THEN
       IF(NSTATM == 0) CALL MSGERR(3,                                    &
              'In stationary mode, the current field also have to be stationary')
       NSTATM = 1
       CALL ININTG('CUR_TBEGINP',IFLBEG(IGRD),'REQ',0)
       CALL ININTG('CUR_DELTINP',IFLINT(IGRD),'REQ',0)
       CALL ININTG('CUR_TENDINP',IFLEND(IGRD),'REQ',0)
       IFLDYN(IGRD) = 1
       IFLTIM(IGRD) = IFLBEG(IGRD)
       IF(IGR2 > 0)THEN
         IFLBEG(IGR2) = IFLBEG(IGRD)                                 
         IFLINT(IGR2) = IFLINT(IGRD)                                 
         IFLEND(IGR2) = IFLEND(IGRD)                                 
         IFLDYN(IGR2) = IFLDYN(IGRD)                                 
         IFLTIM(IGR2) = IFLTIM(IGRD)                                 
       ENDIF
       
       CALL INLOGC('INP_CUR_SERI',INP_CUR_SERI,'REQ','T')
       IF(INP_CUR_SERI)THEN
         INP_CUR_NTIME = (IFLEND(IGRD)-IFLBEG(IGRD)+1)/IFLINT(IGRD)
         ALLOCATE(UXB_GL(MGL,INP_CUR_NTIME)); UXB_GL = 0.0
         ALLOCATE(UYB_GL(MGL,INP_CUR_NTIME)); UYB_GL = 0.0

         ALLOCATE(UXB(MT,INP_CUR_NTIME));  UXB = 0.0
         ALLOCATE(UYB(MT,INP_CUR_NTIME));  UYB = 0.0

	 DO J=1,INP_CUR_NTIME
	   READ(1) (UXB_GL(I,J),UYB_GL(I,J),I=1,MGL)
	 END DO
       ELSE	 
         INP_CUR_NTIME = 1
         ALLOCATE(UXB_GL(MGL,1)); UXB_GL = 0.0
         ALLOCATE(UYB_GL(MGL,1)); UYB_GL = 0.0
	 
         ALLOCATE(UXB(MT,1)); UXB = 0.0
         ALLOCATE(UYB(MT,1)); UYB = 0.0

	 READ(1) (UXB_GL(I,1),UYB_GL(I,1),I=1,MGL)
       END IF	 
     ELSE
       INP_CUR_NTIME = 1
       ALLOCATE(UXB_GL(MGL,1)); UXB_GL = 0.0
       ALLOCATE(UYB_GL(MGL,1)); UYB_GL = 0.0

       ALLOCATE(UXB(MT,1)); UXB = 0.0
       ALLOCATE(UYB(MT,1)); UYB = 0.0

       READ(1) (UXB_GL(I,1),UYB_GL(I,1),I=1,MGL)
     END IF  
     CLOSE(1)   
     
     IF(SERIAL)THEN
       UXB = UXB_GL
       UYB = UYB_GL
     END IF
     
#  if defined (MULTIPROCESSOR)
     IF(PAR)THEN
       DO J=1,INP_CUR_NTIME 
         DO I=1,M
	   UXB(I,J) = UXB_GL(NGID(I),J)
	   UYB(I,J) = UYB_GL(NGID(I),J)
	 END DO
	 DO I=1,NHN
	   UXB(I+M,J) = UXB_GL(HN_LST(I),J)
	   UYB(I+M,J) = UYB_GL(HN_LST(I),J)         
       	 END DO
       END DO
     END IF
#  endif       	 	   
   END IF
   IF(.NOT.ALLOCATED(UXB))    ALLOCATE(UXB(MT,1))
   IF(.NOT.ALLOCATED(UYB))    ALLOCATE(UYB(MT,1))
   
   VARWI = .FALSE.
   INP_WI = .FALSE.
   INP_WI_SERI = .FALSE.

   IF(WIND_ON .AND. WIND_KIND == VRBL) INP_WI = .TRUE.
   IF(INP_WI)THEN
     LWINDR = 2                                                        
     IWIND  = LWINDM                                                   
     IGRD = 5
     IGR1   = 5
     IGR2 = 6
     VARWI  = .TRUE.
     PSNAME = 'WVGRID  '
     INP_WI_SERI = .TRUE.
     
     NSTATM = 1
     CALL SURFACE_WIND2WAVE(IGRD)

     IFLDYN(IGRD) = 1
     IFLTIM(IGRD) = IFLBEG(IGRD)
     IF(IGR2 > 0)THEN
       IFLBEG(IGR2) = IFLBEG(IGRD)                                 
       IFLINT(IGR2) = IFLINT(IGRD)                                 
       IFLEND(IGR2) = IFLEND(IGRD)                                 
       IFLDYN(IGR2) = IFLDYN(IGRD)                                 
       IFLTIM(IGR2) = IFLTIM(IGRD)                                 
     ENDIF
   END IF
   
   CALL INLOGC('INP_FR',INP_FR,'REQ','T')
   IF(INP_FR)THEN
     IGRD = 4
     IGR1   = 4
     VARFR  = .TRUE.
     PSNAME = 'FRICGRID'
     
     CALL INCSTR('INP_FR_NAME',INP_FR_NAME,'REQ',' ')
     INQUIRE(FILE=TRIM(INP_FR_NAME),EXIST=CHECK)
     IF(CHECK)THEN
       OPEN(1,FILE=TRIM(INP_FR_NAME),STATUS='OLD',FORM='UNFORMATTED')
     ELSE
       WRITE(PRINTF,*) TRIM(INP_FR_NAME), ' DOES NOT EXIT.'
       CALL PSTOP
     END IF
    
     CALL INLOGC('INP_FR_STAT',INP_FR_STAT,'REQ','T')
     IF(.NOT. INP_FR_STAT)THEN
       IF(NSTATM == 0) CALL MSGERR(3,                                    &
              'In stationary mode, the friction field also have to be stationary')
       NSTATM = 1
       CALL ININTG('FR_TBEGINP',IFLBEG(IGRD),'REQ',0)
       CALL ININTG('FR_DELTINP',IFLINT(IGRD),'REQ',0)
       CALL ININTG('FR_TENDINP',IFLEND(IGRD),'REQ',0)
       IFLDYN(IGRD) = 1
       IFLTIM(IGRD) = IFLBEG(IGRD)
       IF(IGR2 > 0)THEN
         IFLBEG(IGR2) = IFLBEG(IGRD)                                 
         IFLINT(IGR2) = IFLINT(IGRD)                                 
         IFLEND(IGR2) = IFLEND(IGRD)                                 
         IFLDYN(IGR2) = IFLDYN(IGRD)                                 
         IFLTIM(IGR2) = IFLTIM(IGRD)                                 
       ENDIF
       
       CALL INLOGC('INP_FR_SERI',INP_FR_SERI,'REQ','T')
       IF(INP_FR_SERI)THEN
         INP_FR_NTIME = (IFLEND(IGRD)-IFLBEG(IGRD)+1)/IFLINT(IGRD)
         ALLOCATE(FRIC_GL(MGL,INP_FR_NTIME)); FRIC_GL = 0.0
	 ALLOCATE(FRIC(MT,INP_FR_NTIME));     FRIC    = 0.0
	 DO J=1,INP_FR_NTIME
	   READ(1) (FRIC_GL(I,J),I=1,MGL)
	 END DO
       ELSE	 
         INP_FR_NTIME = 1
         ALLOCATE(FRIC_GL(MGL,1));     FRIC_GL = 0.0
	 ALLOCATE(FRIC(MT,1));         FRIC    = 0.0
	 READ(1) (FRIC_GL(I,1),I=1,MGL)
       END IF	 
     ELSE
       INP_FR_NTIME = 1
       ALLOCATE(FRIC_GL(MGL,1));       FRIC_GL = 0.0
       ALLOCATE(FRIC(MT,1));           FRIC    = 0.0
       READ(1) (FRIC_GL(I,1),I=1,MGL)
     END IF  
     CLOSE(1)
     
     IF(SERIAL)THEN
       FRIC = FRIC_GL
     END IF
     
#    if defined (MULTIPROCESSOR)
     IF(PAR)THEN
       DO J=1,INP_FR_NTIME
         DO I=1,M
	   FRIC(I,J) = FRIC_GL(NGID(I),J)
	 END DO
	 DO I=1,NHN
	   FRIC(I+M,J) = FRIC_GL(HN_LST(I),J)
	 END DO
       END DO
     END IF
#    endif       	            
   END IF
   IF(.NOT.ALLOCATED(FRIC))    ALLOCATE(FRIC(MT,1))
   
   CALL INLOGC('INP_WLEV',INP_WLEV,'REQ','T')
   IF(INP_WLEV)THEN
     IGRD = 7
     IGR1   = 7
     VARWLV = .TRUE.                                                   
     PSNAME = 'WLEVGRID'
     
     CALL INCSTR('INP_WLEV_NAME',INP_WLEV_NAME,'REQ',' ')
     INQUIRE(FILE=TRIM(INP_WLEV_NAME),EXIST=CHECK)
     IF(CHECK)THEN
       OPEN(1,FILE=TRIM(INP_WLEV_NAME),STATUS='OLD',FORM='UNFORMATTED')
     ELSE
       WRITE(PRINTF,*) TRIM(INP_WLEV_NAME), ' DOES NOT EXIT.'
       CALL PSTOP
     END IF
    
     CALL INLOGC('INP_WLEV_STAT',INP_WLEV_STAT,'REQ','T')
     IF(.NOT. INP_WLEV_STAT)THEN
       IF(NSTATM == 0) CALL MSGERR(3,                                    &
              'In stationary mode, the water level also have to be stationary')
       NSTATM = 1
       CALL ININTG('WLEV_TBEGINP',IFLBEG(IGRD),'REQ',0)
       CALL ININTG('WLEV_DELTINP',IFLINT(IGRD),'REQ',0)
       CALL ININTG('WLEV_TENDINP',IFLEND(IGRD),'REQ',0)
       IFLDYN(IGRD) = 1
       IFLTIM(IGRD) = IFLBEG(IGRD)
       IF(IGR2 > 0)THEN
         IFLBEG(IGR2) = IFLBEG(IGRD)                                 
         IFLINT(IGR2) = IFLINT(IGRD)                                 
         IFLEND(IGR2) = IFLEND(IGRD)                                 
         IFLDYN(IGR2) = IFLDYN(IGRD)                                 
         IFLTIM(IGR2) = IFLTIM(IGRD)                                 
       ENDIF
       
       CALL INLOGC('INP_WLEV_SERI',INP_WLEV_SERI,'REQ','T')
       IF(INP_WLEV_SERI)THEN
         INP_WLEV_NTIME = (IFLEND(IGRD)-IFLBEG(IGRD)+1)/IFLINT(IGRD)
         ALLOCATE(WLEVL_GL(MGL,INP_WLEV_NTIME)); WLEVL_GL = 0.0
	 ALLOCATE(WLEVL(MT,INP_WLEV_NTIME));     WLEVL    = 0.0
	 DO J=1,INP_WLEV_NTIME
	   READ(1) (WLEVL_GL(I,J),I=1,MGL)
	 END DO
       ELSE	 
         INP_WLEV_NTIME = 1
         ALLOCATE(WLEVL_GL(MGL,1));    WLEVL_GL = 0.0
	 ALLOCATE(WLEVL(MT,1));        WLEVL    = 0.0
	 READ(1) (WLEVL_GL(I,1),I=1,MGL)
       END IF	 
     ELSE
       INP_WLEV_NTIME = 1
       ALLOCATE(WLEVL_GL(MGL,1));    WLEVL_GL = 0.0
       ALLOCATE(WLEVL(MT,1));        WLEVL    = 0.0
       READ(1) (WLEVL_GL(I,1),I=1,MGL)
     END IF  
     CLOSE(1)
     
     IF(SERIAL)THEN
       WLEVL = WLEVL_GL
     END IF
     
#    if defined (MULTIPROCESSOR)
     IF(PAR)THEN
       DO J=1,INP_WLEV_NTIME
         DO I=1,M
	   WLEVL(I,J) = WLEVL_GL(NGID(I),J)
	 END DO
	 DO I=1,NHN
	   WLEVL(I+M,J) = WLEVL_GL(HN_LST(I),J)
	 END DO
       END DO
     END IF
#    endif       	          
   END IF
   IF(.NOT.ALLOCATED(WLEVL)) ALLOCATE(WLEVL(MT,1))
 
!
!     ------------------------------------------------------------------
!
!     WIND      parameters uniform wind field
!
! ==========================================
!
!   WIND  [vel]  [dir]  [astd]                                            
!
! ==========================================
!
!   CALL INLOGC('WIND',VARWI1,'UNC','T')
    IF(WIND_ON .AND. WIND_KIND == CNSTNT) VARWI1 = .TRUE.
    IF(VARWI1)THEN
     VARWI = .NOT. VARWI1
     LWINDR = 1
     IWIND  = LWINDM
     
     U10  = SQRT(WIND_X*WIND_X+WIND_Y*WIND_Y)
     WDIP = ATAN2(WIND_Y,WIND_X)
     WDIP = WDIP*180.0/PI_W
!
!    *** Convert (if necessary) WDIP from nautical degrees ***         
!    *** to cartesian degrees                              ***         
!
!
     WDIP = DEGCNV (WDIP)
!
     IF(IWIND == 0) IWIND = 4                                         
     ALTMP = WDIP / 360.
     WDIP = PI2_W * (ALTMP - NINT(ALTMP))
     IF(JWX2 <= 1)THEN                                               
       MCMVAR = MCMVAR + 2                                             
       JWX2   = MCMVAR - 1
       JWY2   = MCMVAR
     ENDIF
   ENDIF

!
!  -----------------------------------------------------------------------
!
!  BOUNdary  defining boundary conditions                              
!
!==============================================
!
   NESTING = .FALSE.
   CALL INCSTR('BOUNDARY',BOUNDARY,'REQ',' ')
   
   IF(BOUNDARY /= 'default')THEN

!
!----------------Determine Number of Nodes on Outer Boundary-------------------!
!
     IOBCN_GL_W = IOBCN_GL

     IOBCN_W = 0

     IF(IOBCN_GL_W > 0)THEN

!------------Read in Open Boundary Nodes and Temperature/Salinity Conditions---!

       ALLOCATE(I_OBC_GL_W(IOBCN_GL_W))

       I_OBC_GL_W = I_OBC_GL

!----------------------Make Sure It Is In Global Domain------------------------!

       DO I=1,IOBCN_GL_W
         IF((I_OBC_GL_W(I) > MGL))THEN
           WRITE(IPT,*)'==================ERROR=================================='
           WRITE(IPT,*)'OPEN BOUNDARY NODE NUMBER',I,'IS NOT IN THE'
           WRITE(IPT,*)'GLOBAL DOMAIN'
           WRITE(IPT,*)'CHECK INPUT FILE AND ENSURE OPEN BOUNDARY NODES <= ',MGL
           WRITE(IPT,*)'========================================================='
           CALL PSTOP
         END IF
       END DO

!----------Shift Open Boundary Node List,Type-----------!

       IF(SERIAL)THEN

         IOBCN_W    = IOBCN_GL_W

         ALLOCATE(I_OBC_N_W(IOBCN_W))
         I_OBC_N_W = I_OBC_GL_W
       END IF

#      if defined (MULTIPROCESSOR)
       IF(PAR)THEN
         ALLOCATE(TEMP1(IOBCN_GL_W))

         NCNT = 0
         !!SET UP LOCAL OPEN BOUNDARY NODES
         DO I=1,IOBCN_GL_W
           I1 = NLID( I_OBC_GL_W(I) )
           IF(I1 /= 0)THEN
             NCNT = NCNT + 1
             TEMP1(NCNT) = I1
           END IF
         END DO
         IOBCN_W = NCNT

         IF(NCNT > 0)THEN
           ALLOCATE(I_OBC_N_W(NCNT))
           I_OBC_N_W  = TEMP1(1:NCNT)
         END IF

         DEALLOCATE(TEMP1)
       END IF
#      endif

!adjust ISONB_W according to FVCOM
     ISONB_W = ISONB

     DO I=1,MT
        IF(ISONB(I)/=0) ISONB_W(I)=1 !separate intetior points and outer boundary
     ENDDO

     DO I=1,IOBCN_W
        I1=I_OBC_N_W(I)
        ISONB_W(I1)=2                ! seprate open boundary from outer boundary
     ENDDO
         
     END IF !!IOBCN_GL_W > 0

     CALL SWBOUN ( )  
   ELSE
     IOBCN_W = 0
     ISONB_W = ISONB
     DO I=1,MT
        IF(ISONB_W(I)==2) ISONB_W(I)=1
     ENDDO
           
   END IF
   
!
!  ------------------------------------------------------------------
!
! INITial conditions  Definition of initial conditions for MODE DYNAMIC
!
   IF(KEYWIS('INIT'))THEN                                           
!JQI     CALL INITVA( AC2   , SPCSIG, SPCDIR, KGRPNT, XCGRID, YCGRID, XYTST )
!JQI     IF(STPNOW()) RETURN                                              
!JQI     GO TO 100
   ENDIF
!     ------------------------------------------------------------------
!
!     GEN       the generation of mode                          
!
!     ------------------------------------------------------------------

   CALL ININTG('GEN',GEN,'REQ',0)

!
!     GEN=1        deep water with 1st generation                          
!
! ==========================================
!
! GEN=1     [cf10]  [cf20]  [cf30]  [cf40]  [edmlpm]  [cdrag]  [umin]  [cfpm]
!
! ==========================================
!
   IF(GEN == 1)THEN
     LWINDM = 1                                                        
     IF(LWINDR > 0) IWIND = LWINDM                                 
     IQUAD = 0
     IWCAP = 0

     CALL INREAL('GEN1_CF10',  PWIND(1), 'UNC',0.0)
     CALL INREAL('GEN1_CF20',  PWIND(2), 'UNC',0.0)
     CALL INREAL('GEN1_CF30',  PWIND(3), 'UNC',0.0)
     CALL INREAL('GEN1_CF40',  PWIND(4), 'UNC',0.0)
     CALL INREAL('GEN1_EDMLPM',PWIND(10),'UNC',0.0)
     CALL INREAL('GEN1_CDRAG', PWIND(11),'UNC',0.0)
     CALL INREAL('GEN1_UMIN',  PWIND(12),'UNC',0.0)
     CALL INREAL('GEN1_CFPM',  PWIND(13),'UNC',0.0)
      
!    if [pwtail] is not changed, make it 5
     IF(EQREAL(PWTAIL(1),4.))THEN                                    
       PWTAIL(1) = 5.                                                  
       PWTAIL(3) = PWTAIL(1) + 1.                                      
     ENDIF
   END IF     

!
!     ------------------------------------------------------------------
!
!     GEN2        deep water with 2nd generation                          
!
! ==========================================
!
! GEN2 [cf10] [cf20] [cf30] [cf40] [cf50] [cf60] [edmlpm] [cdrag] [umin] [cfpm]
!
! ==========================================
!
   IF(GEN == 2)THEN
     LWINDM = 2                                                        
     IF(LWINDR > 0) IWIND = LWINDM                                 
     IQUAD = 0
     IWCAP = 0

     CALL INREAL('GEN2_CF10',  PWIND(1), 'UNC',0.0)
     CALL INREAL('GEN2_CF20',  PWIND(2), 'UNC',0.0)
     CALL INREAL('GEN2_CF30',  PWIND(3), 'UNC',0.0)
     CALL INREAL('GEN2_CF40',  PWIND(4), 'UNC',0.0)
     CALL INREAL('GEN2_CF50',  PWIND(5), 'UNC',0.0)
     CALL INREAL('GEN2_CF60',  PWIND(6), 'UNC',0.0)
     CALL INREAL('GEN2_EDMLPM',PWIND(10),'UNC',0.0)
     CALL INREAL('GEN2_CDRAG', PWIND(11),'UNC',0.0)
     CALL INREAL('GEN2_UMIN',  PWIND(12),'UNC',0.0)
     CALL INREAL('GEN2_CFPM',  PWIND(13),'UNC',0.0)
      
!    if [pwtail] is not changed, make it 5
     IF(EQREAL(PWTAIL(1),4.))THEN                                    
       PWTAIL(1) = 5.                                                  
       PWTAIL(3) = PWTAIL(1) + 1.                                      
     ENDIF
   END IF     

!
!     ------------------------------------------------------------------
!
!     GEN3        deep water with 3d generation                           
!
! ======================================================================================
!
!       |   JANSsen [cds1] [delta]                         |              
! GEN3 <                                                    >     &
!       | ->KOMen   [cds2] [stpm] [powst] [delta] [powk]   |              
!       |                                                  |
!       |   YAN     (NOT documented)                       |
!       |                                                  |
!       |   WESTHuysen [cds2] [br] [p0] [powst] [powk]     |              
!
!  (QUADrupl [iquad] [limiter] [lambda] [cnl4] [csh1] [csh2] [csh3]) (AGROW [a])
!
! ======================================================================================
!
   IF(GEN == 3)THEN
     GROWTH = 'KOM'
     CALL INCSTR('GROWTH',GROWTH,'STA','KOM')
     
     IF(GROWTH == 'JANS')THEN
       LWINDM = 4
       IF(LWINDR > 0) IWIND = LWINDM
!      *** whitecapping according to Janssen (1989, 1991) ***
       IWCAP = 2
       CALL INREAL('GEN3_JANS_CDS1', PWCAP(3),'UNC',0.0)
       CALL INREAL('GEN3_JANS_DELTA',PWCAP(4),'UNC',0.0)
       
!
!      Recalculate coefficientes that are actually used in the         
!      whitecapping routines                                           
       PWCAP(1)  = PWCAP(3) * (PWCAP(2) ** PWCAP(9))                   
       PWCAP(10) = PWCAP(4)                                            
       JUSTAR = MCMVAR+1                                               
       JZEL   = MCMVAR+2                                               
       JCDRAG = MCMVAR+3                                               
       JTAUW  = MCMVAR+4                                               
       MCMVAR = MCMVAR+4                                               
!      if [pwtail] is not changed, make it 5
       IF(EQREAL(PWTAIL(1),4.))THEN                                  
         PWTAIL(1) = 5.                                                
         PWTAIL(3) = PWTAIL(1) + 1.                                    
       ENDIF
   
     ELSE IF(GROWTH == 'YAN')THEN
!      option not documented in user manual
       LWINDM = 5                                                      
       IF(LWINDR > 0) IWIND = LWINDM                               
     ELSE IF(GROWTH == 'WESTH')THEN                                    
!      *** wind according to (adapted) Yan (1987) ***                  
       LWINDM = 5                                                      
       IF(LWINDR > 0) IWIND = LWINDM                               
!      *** whitecapping according to Alves and Banner (2003) ***       
       IWCAP = 7                                                       

       CALL INREAL('GEN_WESTH_CDS2', PWCAP(1), 'STA',5.0E-5)
       CALL INREAL('GEN_WESTH_BR',   PWCAP(12),'STA',1.75E-3)                
       CALL INREAL('GEN_WESTH_P0',   PWCAP(10),'STA',4.)                     
       CALL INREAL('GEN_WESTH_POWST',PWCAP(9), 'STA',0.)                     
       CALL INREAL('GEN_WESTH_POWK', PWCAP(11),'STA',0.)                     
     ELSE IF(GROWTH == 'KOM')THEN
       LWINDM = 3                                                      
       IF(LWINDR > 0) IWIND = LWINDM                               
!      *** whitecapping according to Komen et al. (1984) ***
       IWCAP = 1

       CALL INREAL('GEN3_KOM_CDS2',PWCAP(1),'UNC',0.0)
       CALL INREAL('GEN3_KOM_STPM',PWCAP(2),'UNC',0.0)
     ELSE
       WRITE(*,*)
       CALL PSTOP      
     ENDIF

!    ***  parameters nonlinear 4 wave interactions ***
     QUAD = .TRUE.
     CALL INLOGC('QUAD',QUAD,'STA','T')
     
     IF(QUAD)THEN
       CALL ININTG('IQUAD',IQUAD,'UNC',0)
      
!      *** if quadruplets are activated then LIMITER = 0.1 ***
!      *** the standard value, however, is 1.e20           ***

       CALL INREAL('LAMBDA',PQUAD(1),'UNC',0.0)
       CALL INREAL('CNL4',  PQUAD(2),'UNC',0.0)
       CALL INREAL('CSH1',  PQUAD(3),'UNC',0.0)
       CALL INREAL('CSH2',  PQUAD(4),'UNC',0.0)
       CALL INREAL('CSH3',  PQUAD(5),'UNC',0.0)
     ELSE
       IQUAD = 0
     END IF 
     
     AGROW = .FALSE.
     CALL INLOGC('AGROW',AGROW,'STA','F')

     IF(AGROW)THEN
       CALL INREAL('A',PWIND(31),'STA',0.0015)
     ELSE
       IF(NSTATM == 1 .AND. ICOND == 0) ICOND = 1
     END IF  
   END IF
!
!     ------------------------------------------------------------------
!
!     WCAP        parameters whitecapping
!
! =============================================================
!
!        | ->KOMen    [cds2] [stpm] [powst] [delta] [powk]  |             
!        |                                                  |
!        |   JANSsen  [cds1]  [delta] [pwtail]              |
! WCAP  <                                                    >    (NOT documented)
!        |   LHIG     [cflhig]                              |
!        |                                                  |
!        |   BJ       [bjstp] [bjalf]                       |
!        |                                                  |
!        |   KBJ      [bjstp] [bjalf] [kconv]               |
!        |                                                  |             
!        |   CSM      [cst] [pow]                           |             
!        |                                                  |
!        |   AB       [cds2] [br] [p0] [powst] [powk]       |             
!
! =============================================================
!
   WCAP = 'KOM '
   CALL INCSTR('WCAP',WCAP,'STA','KOM')

   IF(WCAP == 'KOM')THEN
!    *** whitecapping according to Komen et al. (1984) ***
     IWCAP = 1
     CALL INREAL('WCAP_KOM_CDS2', PWCAP(1), 'UNC',0.0)
     CALL INREAL('WCAP_KOM_STPM', PWCAP(2), 'UNC',0.0)
     CALL INREAL('WCAP_KOM_POWST',PWCAP(9), 'UNC',0.0)
     CALL INREAL('WCAP_KOM_DELTA',PWCAP(10),'UNC',0.0)
     CALL INREAL('WCAP_KOM_POWK', PWCAP(11),'UNC',0.0)
   ELSE IF(WCAP == 'JANS')THEN
!    *** whitecapping according to Janssen (1989, 1991) ***
     IWCAP = 2
     CALL INREAL('WCAP_JANS_CDS1', PWCAP(3),'UNC',0.0)  
     CALL INREAL('WCAP_JANS_DELTA',PWCAP(4),'UNC',0.0)  
!
!    Recalculate coefficients that are actually used in the          
!    whitecapping routines                                           
     PWCAP(1)  = PWCAP(3) * (PWCAP(2) ** PWCAP(9))                   
     PWCAP(10) = PWCAP(4)                                            
     JUSTAR = MCMVAR+1                                               
     JZEL   = MCMVAR+2                                               
     JCDRAG = MCMVAR+3                                               
     JTAUW  = MCMVAR+4                                               
     MCMVAR = MCMVAR+4                                               

     CALL INREAL('WCAP_JANS_PWTAIL',PWTAIL(1),'STA',5.0)
     IF(PWTAIL(1) <= 1.) CALL MSGERR (3, 'Incorrect PWTAIL')        
     PWTAIL(3) = PWTAIL(1) + 1.                                      
   ELSE IF(WCAP == 'LHIG')THEN
!    *** whitecapping according to Longuett Higgins ***
     IWCAP = 3
     PWCAP(5) = 1.0
     CALL INREAL('WCAP_LHIG_CFLHIG',PWCAP(5),'UNC',0.0)
   ELSE IF(WCAP == 'BJ')THEN
!    *** whitecapping according to Battjes/Janssen formulation ***
     IWCAP = 4
     PWCAP(6) = 0.88
     CALL INREAL('WCAP_BJ_BJSTP',PWCAP(6),'UNC',0.0)
     PWCAP(7) = 1.0
     CALL INREAL('WCAP_BJ_BJALF',PWCAP(7),'UNC',0.0)
   ELSE IF(WCAP == 'KBJ')THEN
!    *** whitecapping according to a combination of Komen et al ***
!    *** and Battjes/Janssen                                    ***
     IWCAP = 5
     PWCAP(6) = 0.88
     CALL INREAL('WCAP_KBJ_BJSTP',PWCAP(6),'UNC',0.0)
     PWCAP(7) = 1.0
     CALL INREAL('WCAP_KBJ_BJALF',PWCAP(7),'UNC',0.0)
     PWCAP(8) = 0.75
     CALL INREAL('WCAP_KBJ_KCONV',PWCAP(8),'UNC',0.0)
   ELSE IF(WCAP == 'CSM')THEN                                    
!    *** whitecapping according to cumulative steepness method  ***  
     IWCAP = 6                                                       
     PWCAP(12) = 4.0
     CALL INREAL('WCAP_CSM_CST',PWCAP(12),'UNC',0.0)
     PWCAP(13) = 2.0
     CALL INREAL('WCAP_CSM_POW',PWCAP(13),'UNC',0.0)
   ELSE IF(WCAP == 'AB')THEN                                     
!    *** whitecapping according to Alves and Banner (2003) ***       
     IWCAP = 7                                                       
     CALL INREAL('WCAP_AB_CDS2',PWCAP(1), 'STA',5.0E-5)
     CALL INREAL('WCAP_AB_BR',  PWCAP(12),'STA',1.75E-3)
     CALL INREAL('WCAP_AB_PO',  PWCAP(10),'STA',4.0)
     CALL INREAL('WCAP_AB_POWST',PWCAP(9),'STA',0.0)
     CALL INREAL('WCAP_AB_POWK',PWCAP(11),'STA',0.0)
   ELSE IF(WCAP == 'OFF')THEN
     IWCAP = 0
   ELSE
!JQI          CALL WRNKEY
   ENDIF

!     ------------------------------------------------------------------
!     ***  parameters nonlinear 4 wave interactions ***
     QUAD = .TRUE.
     CALL INLOGC('QUAD',QUAD,'STA','T')
     
     IF(QUAD)THEN
!
! ===================================================================
!
!  QUADrupl [iquad] [limiter] [lambda] [cnl4] [csh1] [csh2] [csh3]        
!
! ===================================================================
!
       CALL ININTG('IQUAD',IQUAD,'UNC',0)
      
!      *** if quadruplets are activated then LIMITER = 0.1 ***
!      *** the standard value, however, is 1.e20           ***
       CALL INREAL('LIMITER',PNUMS(20),'STA',0.1)           
       CALL INREAL('LAMBDA', PQUAD(1), 'UNC',0.0)
       CALL INREAL('CNL4',   PQUAD(2), 'UNC',0.0)
       CALL INREAL('CSH1',   PQUAD(3), 'UNC',0.0)
       CALL INREAL('CSH2',   PQUAD(4), 'UNC',0.0)
       CALL INREAL('CSH3',   PQUAD(5), 'UNC',0.0)
     ELSE
       IQUAD = 0
     ENDIF

! ----------------------------------------------------------------    
!                            | CNL4 < [cnl4] >               |         
!  MDIA LAMbda < [lambda] > <                                 >        
!                            | CNL4_12 < [cnl4_1] [cnl4_2] > |         
! ----------------------------------------------------------------    

   MDIA = .FALSE.
   CALL INLOGC('MDIA',MDIAL,'UNC','F')

   IF(MDIAL)THEN                                            
     IF(ALLOCATED(LAMBDA))THEN                                       
       DEALLOCATE (LAMBDA,CNL4_1,CNL4_2)                               
     ENDIF                                                             
     ALLOCATE (RLAMBDA(1000))                                          

     LAM = .FALSE.
     CALL INLOGC('MDIA_LAM',LAM,'STA','F')
     
     IF(LAM)THEN                                           
       MORE = .TRUE.                                                   
       ILAMBDA = 0                                                     
       DO WHILE (MORE)                                                 
         CALL INREAL('MDIA_LAMBDA',RLAMBDA(ILAMBDA+1),'STA',-1.)            
         IF(RLAMBDA(ILAMBDA+1) < 0.) THEN                            
           MORE = .FALSE.                                              
         ELSE                                                          
           ILAMBDA = ILAMBDA + 1                                       
         ENDIF                                                         
       ENDDO                                                           
       MDIA = ILAMBDA                                                  
       ALLOCATE (LAMBDA(MDIA))                                         
       LAMBDA(1:MDIA) = RLAMBDA(1:MDIA)                                
       DEALLOCATE (RLAMBDA)                                            
     ELSE                                                              
!JQI   CALL WRNKEY                                                     
     END IF                                                            
     ALLOCATE (CNL4_1(MDIA),CNL4_2(MDIA))                              
     CALL INCSTR('MDIA_CNL4C',CNL4,'UNC',' ')
     IF(CNL4 == 'CNL4_12')THEN                                       
       DO ICNL4=1,MDIA
         CALL INREAL('MDIA_CNL4_1',CNL4_1(ICNL4),'REQ',0.)                  
         CALL INREAL('MDIA_CNL4_2',CNL4_2(ICNL4),'REQ',0.)                  
       END DO                                                          
     ELSE IF(CNL4 == 'CNL4')THEN                                      
       DO ICNL4=1,MDIA                                                 
         CALL INREAL('MDIA_CNL4',CNL4_1(ICNL4),'REQ',0.)                    
       END DO                                                          
       CNL4_2 = CNL4_1                                                 
     ELSE                                                              
!JQI          CALL WRNKEY                                                     
     END IF                                                            
     CNL4_1 = CNL4_1 * ((2.*PI_W)**9)                                    
     CNL4_2 = CNL4_2 * ((2.*PI_W)**9)                                    
   ENDIF                                                               

!
! ------------------------------------------------------------------
!
!     BREAK     parameters surf breaking
!
! ============================================================
!
!           | -> CONstant [alpha] [gamma]                                      |
! BREaking <                                                                    >
!           |    VARiable [alpha] [gammin] [gammax] [gamneg] [coeff1] [coeff2] |
!
! ============================================================
!
  BRE = .FALSE.
  CALL INLOGC('BRE',BRE,'STA','F')
  IF(BRE)THEN
    CALL INCSTR('CONSTANT',CONSTANT,'STA','CON')

    IF(CONSTANT == 'CON')THEN
      ISURF = 1
      CALL INREAL('BRE_CON_ALPHA',PSURF(1),'STA',1.0)
      CALL INREAL('BRE_CON_GAMMA',PSURF(2),'STA',0.73)
    ELSE IF(CONSTANT == 'VAR')THEN                                      
      ISURF = 2
      CALL INREAL('BRE_VAR_ALPHA', PSURF(1),'STA',1.5)
      CALL INREAL('BRE_VAR_GAMMIN',PSURF(4),'STA',0.55)
      CALL INREAL('BRE_VAR_GAMMAX',PSURF(5),'STA',0.81)
      CALL INREAL('BRE_VAR_GAMNEG',PSURF(6),'STA',0.73)
      CALL INREAL('BRE_VAR_COEFF1',PSURF(7),'STA',0.88)
      CALL INREAL('BRE_VAR_COEFF2',PSURF(8),'STA',0.012)
    ENDIF
  ELSE
    ISURF = 0
  ENDIF
!
! ------------------------------------------------------------------
!
!     FRICTION  setting parameters for bottom friction
!
! ===============================================
!
!                   | -> JONSWAP       [cfjon]
!                   |
!   FRICTION       <     COLLINS       [cfw]    &
!                   |
!                   |                  [cfc]    (NOT documented)
!                   |
!                   |    MADSEN        [kn]
!
! ===============================================
!
   FRICL = .FALSE.
   CALL INLOGC('FRICTION',FRICL,'UNC','F')
   
   IF(FRICL)THEN
     IBOT = 1
     CALL INCSTR('FRIC_FORM',FRIC_FORM,'STA','JON')
     IF(FRIC_FORM == 'JON')THEN
       IBOT = 1
       CALL INREAL('CFJON',PBOT(3),'UNC',0.)
     ELSE IF(FRIC_FORM == 'COLL')THEN                                     
       IBOT = 2
       CALL INREAL('CFW',PBOT(2),'UNC',0.)
       CALL INREAL('CFC',PBOT(1),'UNC',0.)
     ELSE IF(FRIC_FORM == 'MAD')THEN
       IBOT = 3
       CALL INREAL('KN',PBOT(5),'UNC',0.)
     ELSE
!JQI       CALL WRNKEY
     ENDIF
   ENDIF

!
! ------------------------------------------------------------------
!
!     ***  parameters nonlinear 3 wave interactions
!
! ============================================================
!
!  TRIad [trfac] [cutfr] [urcrit] [urslim]                                
!
! ============================================================
!
   TRI = .FALSE.
   CALL INLOGC('TRIAD',TRI,'UNC','F')
   
   IF(TRI)THEN
     ITRIAD = 1                                                        
     CALL INREAL('TRFAC' ,PTRIAD(1),'UNC',0.0)
     CALL INREAL('CUTFR', PTRIAD(2),'STA',2.5)                      
     CALL INREAL('URCRIT',PTRIAD(4),'UNC',0.0)                     
     CALL INREAL('URSLIM',PTRIAD(5),'STA',0.01)                    
!    reserve space in array Compda for Ursell number                   
     IF(JURSEL <= 1)THEN
       MCMVAR = MCMVAR + 1
       JURSEL = MCMVAR                                                 
     ENDIF
   ENDIF

!
! ------------------------------------------------------------------
!
!     LIM     setting parameters in conjunction with action limiter
!
! ============================================================
!
!   LIMiter [ursell] [qb]                                                 
!
! ============================================================
!
   LIM = .FALSE.
   CALL INLOGC('LIMITER',LIM,'UNC','F')
   
   IF(LIM)THEN                                            
     CALL INREAL('URSELL',PTRIAD(3),'UNC',0.)                          
     CALL INREAL('QB'    ,PNUMS(28),'UNC',0.)                          
   ENDIF                                                               

!
! ------------------------------------------------------------------
!
!     *** OBSTACLE   Definition of obstacles in comp grid. ***            
!
! ============================================================
!
!             |  TRANSm [trcoef]              |
! OBSTacle   <                                 >             &
!             |  DAM    [hgt] [alpha] [beta]  |
!                                                                         
!                        | -> RSPEC                    |                  
!      ( REFLec [reflc] <                              > )   &
!                        |    RDIFF [POWS][POWD][Kdif] |                  
!
!             LINe < [xp] [yp] >                                          
!
! ============================================================
!
   OBST = .FALSE.
   CALL INLOGC('OBSTACLE',OBST,'UNC','F')
   
   IF(OBST)THEN
!
     ALLOCATE(OBSTMP)                                                  
     OBSTMP%TRCOEF(1) = 0.                                             
     OBSTMP%TRCOEF(2) = 0.                                             
     OBSTMP%TRCOEF(3) = 0.                                             
     OBSTMP%RFCOEF(1) = 0.                                             
     OBSTMP%RFCOEF(2) = 0.                                             
     OBSTMP%RFCOEF(3) = 0.                                             
     OBSTMP%RFCOEF(4) = 0.                                             
     OBSTMP%RFCOEF(5) = 0.                                             
     OBSTMP%RFCOEF(6) = 0.                                             
     OBSTMP%RFCOEF(7) = 0.                                             
     OBSTMP%RFCOEF(8) = 0.                                             
!
!    data concerning transmission of energy over/through the obstacle
!
     CALL INCSTR('OBST_TYPE',OBST_TYPE,'UDF',' ')

     IF(OBST_TYPE == 'TRANS')THEN
       ITRAS = 0
       CALL INREAL('TRCOEF',TRCF,'UDF', 0.)
       IF((TRCF < 0.) .OR. (TRCF > 1.))THEN                        
         CALL MSGERR(3,'Transmission coeff. [trans] is not allowed ')  
         CALL MSGERR(3,'to be greater than 1 or smaller than 0!  ')    
       ENDIF                                                           
       OBSTMP%TRCOEF(1) = TRCF                                         
     ELSE IF(OBST_TYPE == 'DAM')THEN
       ITRAS = 1
       CALL INREAL('HGT'  , HGT , 'UDF', 0.  )                        
       CALL INREAL('ALPHA', OGAM, 'STA', 2.6 )                        
       CALL INREAL('BETA' , OBET, 'STA', 0.15)                        
       OBSTMP%TRCOEF(1) = HGT                                          
       OBSTMP%TRCOEF(2) = OGAM                                         
       OBSTMP%TRCOEF(3) = OBET                                         
     ELSE
!      if no transmission options are activated, there will be 0 transmission
       ITRAS = 0
       TRCF  = 0.                                                      
       OBSTMP%TRCOEF(1) = TRCF                                         
     ENDIF
     OBSTMP%TRTYPE = ITRAS                                             
!                                                                      
!    data on reflection by the obstacle, activated by keyword REFL     
!
     REFL = .FALSE.
     CALL INLOGC('REFLEC',REFL,'DFL','T')

     IF(REFL)THEN                                         
       LREF = 1                                                     
       IF(.NOT.FULCIR)THEN                                           
         CALL MSGERR(3,'Reflections will only be calculated if     ')  
         CALL MSGERR(3,'the spectral directions cover the full     ')  
         CALL MSGERR(3,'circle.                                    ')  
       ENDIF                                                           
       CALL INREAL('REFLC',REF0,'STA',1.)                          
       DUM = REF0*REF0+TRCF*TRCF                                       
       IF(ITRAS == 0 .AND. ((DUM < 0.) .OR. (DUM > 1.)))THEN         
         CALL MSGERR(3,'Kt^2 + Kr^2 > 1 or Kt^2 + Kr^2 < 0! ')        
       ENDIF                                                          
       OBSTMP%RFCOEF(1) = REF0                                        
!                                                                     
       CALL INCSTR('RDIFF',RDIFF,'UDF',' ')

       IF(RDIFF == 'RDIFF')THEN                                       
         LREFDIFF = 1                                                
         CALL INREAL('POWS',POWS,'UDF',1.)                       
         DUM = MOD(POWS,1.)                                          
         IF(POWS < 0.)THEN                                           
           CALL MSGERR(3,'Power POWS is not a positive number! ')    
         ENDIF                                                       
         IF(DUM /= 0.)THEN                                         
           CALL MSGERR(3,'Power POWS is not an integer number! ')     
         ENDIF                                                         
         OBSTMP%RFCOEF(2) = POWS                                       
!                                                                      
         CALL INREAL('POWD',POWD,'STA',1.)                         
         DUM = MOD(POWD,1.)                                            
         IF(POWD < 0.)THEN                                           
           CALL MSGERR(3,'Power POWD is not a positive number! ')     
         ENDIF                                                        
         IF(DUM /= 0.)THEN                                         
           CALL MSGERR(3,'Power POWD is not an integer number! ')      
         ENDIF                                                         
         OBSTMP%RFCOEF(3) = POWD                                       
         CALL INREAL ('KDIF', KDIF, 'STA', 0.)                         
         IF((KDIF < 0.) .OR. (KDIF > 1.))THEN                        
           CALL MSGERR(3,'KDIF > 1 or KDIF < 0! ')                     
         ENDIF                                                         
         OBSTMP%RFCOEF(4) = KDIF                                       
       ELSE IF(RDIFF == 'RSPEC')THEN                                                           
         LREFDIFF = 0                                                  
       ENDIF                                                           
       OBSTMP%RFTYP2 = LREFDIFF                                        
!
       RFD = .FALSE.
       CALL INLOGC('RFD',RFD,'DFL','T')

       IF(RFD)THEN                                         
         LRFRD = 1                                                     
         CALL INREAL('FD1',FD1,'REQ', 0.565)                        
         CALL INREAL('FD2', FD2, 'REQ', 0.75)                         
         CALL INREAL('FD3', FD3, 'REQ', -21.)                         
         CALL INREAL('FD4', FD4, 'REQ', 0.066)                        
         IF(FD3 >= 0.)THEN                                         
           CALL MSGERR(3,'Positive frequency dependency! ')            
         ENDIF                                                         
         OBSTMP%RFCOEF(5) = FD1                                        
         OBSTMP%RFCOEF(6) = FD2                                        
         OBSTMP%RFCOEF(7) = FD3                                        
         OBSTMP%RFCOEF(8) = FD4                                        
       ELSE                                                            
         LRFRD = 0                                                     
       ENDIF                                                           
       OBSTMP%RFTYP3 = LRFRD                                           

       IF((REF0 < 0.) .OR. (REF0 > 1.))THEN                        
         CALL MSGERR(3,'Reflection coeff. [reflc] is not allowed ')    
         CALL MSGERR(3,'to be greater than 1 or smaller than 0!  ')    
       ENDIF                                                           
     ELSE
!      if there is no keyword REFL, there will be no reflection        
       LREF = 0
     ENDIF                                                             
     OBSTMP%RFTYP1 = LREF                                              
!
!    location of obstacles
!
!    *** NUMCOR : Number of corners ***
     NUMCOR = 0
     CALL INLOGC('LINE',LINE,'REQ','T')
     IF(LINE)THEN
       FRST%X = 0.                                                     
       FRST%Y = 0.                                                     
       NULLIFY(FRST%NEXTXY)                                            
       CURR => FRST                                                    
       DO
!        read coordinates of one corner point of the obstacle
!JQI         CALL READXY ('XP', 'YP', XP, YP, 'REP', -1.E10, -1.E10)
!JQI         IF(XP < -.9E10) GOTO 101
!JQI         NUMCOR = NUMCOR + 1
!JQI         ALLOCATE(TMP)                                                 
!JQI         TMP%X = XP                                                    
!JQI         TMP%Y = YP                                                    
!JQI         NULLIFY(TMP%NEXTXY)                                           
!JQI         CURR%NEXTXY => TMP                                            
!JQI         CURR => TMP                                                   
       ENDDO
     ENDIF
101  OBSTMP%NCRPTS = NUMCOR                                            
!    store coordinates for corner points in array of obstacle data     
     ALLOCATE(OBSTMP%XCRP(NUMCOR))                                     
     ALLOCATE(OBSTMP%YCRP(NUMCOR))                                     
     CURR => FRST%NEXTXY                                               
     DO JJ = 1, NUMCOR                                                 
       OBSTMP%XCRP(JJ) = CURR%X                                      
       OBSTMP%YCRP(JJ) = CURR%Y                                      
       CURR => CURR%NEXTXY                                           
     END DO                                                            
     DEALLOCATE(TMP)                                                   
     NULLIFY(OBSTMP%NEXTOBST)                                          
     IF(NUMCOR <= 1)THEN
!JQI       CALL MSGERR(1,'No corner points for obstacle were found')
     ELSE
!      *** NUMOBS : Number of obstacles ***
       NUMOBS = NUMOBS + 1
       IF(.NOT.LOBST)THEN                                          
         FOBSTAC = OBSTMP                                             
         COBST => FOBSTAC                                             
         LOBST = .TRUE.                                               
       ELSE                                                           
         COBST%NEXTOBST => OBSTMP                                     
         COBST => OBSTMP                                              
       END IF                                                         
     ENDIF
!
   ENDIF

!
!  ------------------------------------------------------------------
!  SETUP     include wave-induced set-up in SWAN calculation
!
! ============================================================
!
!   SETUP [supcor]                                                        
!
! ============================================================
!
   SETUP = .FALSE.
   LSETUP = 0                                                        
   CALL INLOGC('SETUP',SETUP,'UNC','F')
   
   IF(SETUP)THEN  
     IF(KSPHER /= 0)                          &
       CALL MSGERR (2,'setup will not be compute correctly with spherical coordinates')  
     LSETUP = 1                                                        
     CALL INREAL ('SUPCOR',PSETUP(2),'UNC',0.)                      
!    *** set pointers for setup and saved depth in array COMPDA ***    
     JSETUP = MCMVAR + 1                                               
     JDPSAV = MCMVAR + 2                                               
     MCMVAR = MCMVAR + 2                                               
   ENDIF  

!
!  ------------------------------------------------------------------  
!
!  DIFFRac   include diffraction approximation                         
!                                                                         
! =============================================================           
!                                                                         
!   DIFFRac  [idiffr]  [smpar]  [smnum]  [cgmod]                          
!                                                                         
! =============================================================           
!                                                                         
   DIFFR = .FALSE.
   CALL INLOGC('DIFFRAC',DIFFR,'UNC','F')
   
   IF(DIFFR)THEN                                          
     CALL ININTG('IDIFFR', IDIFFR, 'STA', 1)                          
     CALL INREAL('SMPAR', PDIFFR(1), 'STA', 0.0)                      
     CALL INREAL('SMNUM', PDIFFR(2), 'STA', 0.)                       
     CALL INREAL('CGMOD', PDIFFR(3), 'STA', 1.)                       
   ENDIF                                                               

!
!  ------------------------------------------------------------------
!
!     OFF      switching standard options off
!
! ============================================================
!
!        |  BREaking
!        |
!        |  WCAPping
!        |
!        |  REFrac
! OFF   <
!        |  FSHift
!        |
!        |  QUADrupl
!        |
!        |  WINDGrowth                                                    
!        |
!        |  BNDCHK                                                        
!
! ============================================================
!
   OFF = .FALSE.
   CALL INLOGC('OFF',OFF,'UNC','F')
   
   IF(OFF)THEN
     CALL INCSTR('OFF_REF',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') IREFR = 0
     CALL INCSTR('OFF_FSH',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') ITFRE = 0
     CALL INCSTR('OFF_BRE',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') ISURF = 0
     CALL INCSTR('OFF_WCAP',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') IWCAP = 0
     CALL INCSTR('OFF_QUAD',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') IQUAD = 0
     CALL INCSTR('OFF_WINDG',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') IWIND = 0
     CALL INCSTR('OFF_BNDCHK',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') BNDCHK = .FALSE.        !switch off checking of Hs on boundary  
     CALL INCSTR('OFF_RESCALE',OFF_TYPE,'STA','YES')
     IF(OFF_TYPE == 'YES') BRESCL = .FALSE.        !switch off rescaling solution
   END IF
!
!  ------------------------------------------------------------------
! =======================================================================
!
!            |    BSBT                      |
!   PROP    <                     | SEC |   |
!            |    GSE  [waveage] <  MIN  >  |
!            |                    |  HR |   |
!            |                    | DAY |   |
!
!            FLUXLIM
!
! =======================================================================
!
   PROP = .FALSE.
   CALL INLOGC('PROP',PROP,'UNC','F')
   IF(PROP)THEN                                           
     CALL INCSTR('PROP_CHOICE', PROP_CHOICE,'UDF',' ')

     IF(PROP_CHOICE == 'BSBT')THEN                                         
       PROPSN = 1                                                    
       PROPSS = 1                                                    
     ELSE IF(PROP_CHOICE == 'GSE')THEN                                     
       IF(PROPSN /= 3)THEN                                           
         CALL MSGERR(2,'Anti-GSE only allowed for S&L scheme.')    
       ENDIF                                                         
       CALL ININTV('WAVEAGE', WAVAGE, 'STA', 0.)                     
     ENDIF                                                           

     FLUXLIM = .FALSE.
     CALL INLOGC('FLUXLIM',FLUXLIM,'DFL','T')
     
     IF(FLUXLIM)THEN                                      
       PROPFL   = 1                                                  
       PNUMS(6) = 0.                                                 
     END IF                                                          
   ENDIF                                                             

!
!     ------------------------------------------------------------------
!
!     NUM       setting numerical parameters for schemes, solvers, etc.
!
! =====================================================================
!
!             | -> ACCUR [drel] [dhoval] [dtoval] [npnts] |               
!   NUMeric (<                                             >          &   
!             |    STOPC [dabs] [drel] [curvat] [npnts]   |               
!
!                    | -> STAT  [mxitst] [alfa] |                         
!                   <                            >  [limiter]   )     &   
!                    | NONSTat  [mxitns]        |
!
!           ( DIRimpl [cdd] [cdlim]                                )  &
!
!
!           | -> SIGIMpl [css] [eps2] [outp] [niter]               |
!          (<                                                      >) &
!           |    SIGEXpl [cfl]                 (NOT documented)    |
!           |                                                      |
!           |    FIL     [diffc]               (NOT documented)    |
!
!
!           ( SETUP [eps2] [outp] [niter]                          )      
!
! =====================================================================
!
   NUM = .FALSE.
   CALL INLOGC('NUMERIC',NUM,'UNC','F')
   
   IF(NUM)THEN
     CALL INCSTR('NUM_CHOICE',NUM_CHOICE,'UDF',' ')

!    *** accuracy and criterion to terminate the iteration ***
     IF(NUM_CHOICE == 'ACCUR')THEN
       PNUMS(21) = 0.                                                  
       CALL INREAL ('DREL'   , PNUMS(1) , 'UNC', 0.)
       CALL INREAL ('DHOVAL' , PNUMS(15), 'UNC', 0.)                   
       CALL INREAL ('DTOVAL' , PNUMS(16), 'UNC', 0.)                   
       CALL INREAL ('NPNTS'  , PNUMS(4) , 'UNC', 0.)

       CALL INCSTR('STAT',STATL,'UDF',' ')

       IF(STATL == 'STAT')THEN                                       
         CALL ININTG ('MXITST' , MXITST   , 'REQ', 0 )                 
         CALL INREAL ('ALFA'   , PNUMS(30), 'UNC', 0.)                 
       ELSE IF(STATL == 'ITERMX')THEN                                
         CALL ININTG ('MXITST' , MXITST   , 'REQ', 0 )                 
         MXITNS = MXITST                                               
         CALL MSGERR (1, '[itermx] is replaced; see user manual')      
       ELSE IF(STATL == 'NONST')THEN                                 
         CALL ININTG ('MXITNS' , MXITNS   , 'REQ', 0 )                 
       ENDIF
       CALL INREAL ('LIMITER', PNUMS(20), 'UNC', 0.)                   
     ELSE IF(NUM_CHOICE == 'STOPC')THEN                                   
       PNUMS(21) = 1.                                                  
       CALL INREAL ('DABS'   , PNUMS(2) , 'STA', 0.00)                 
       CALL INREAL ('DREL'   , PNUMS(1) , 'STA', 0.01)                 
       CALL INREAL ('CURVAT' , PNUMS(15), 'STA', 0.005)                
       CALL INREAL ('NPNTS'  , PNUMS(4) , 'UNC', 0.)                   

       CALL INCSTR('STAT',STATL,'UDF',' ')

       IF(STATL == 'STAT')THEN                                       
         CALL ININTG ('MXITST' , MXITST   , 'STA', 50)                 
         CALL INREAL ('ALFA'   , PNUMS(30), 'UNC', 0.)                 
       ELSE IF(STATL == 'ITERMX')THEN                                
         CALL ININTG ('MXITST' , MXITST   , 'REQ', 0 )                 
         MXITNS = MXITST                                               
         CALL MSGERR (1, '[itermx] is replaced; see user manual')      
       ELSE IF(STATL == 'NONST')THEN                                 
         CALL ININTG ('MXITNS' , MXITNS   , 'REQ', 0 )                 
       ENDIF                                                           
       CALL INREAL ('LIMITER', PNUMS(20), 'UNC', 0.)                   
     END IF
!
!    *** numerical scheme in directional space (standard  ***
!    *** option implicit scheme)                          ***
!
     DIR = .FALSE.
     CALL INLOGC('DIRIMPL',DIR,'DFL','T')

     IF(DIR)THEN
       CALL INREAL ('CDD'    , PNUMS(6) , 'UNC', 0.)
       CALL INREAL ('CDLIM'  , PNUMS(17), 'UNC', 0.)                   
       IF(PNUMS(17) < 0.) IREFR = 1                                  
       IF(PNUMS(17) > 0.) IREFR = -1                                 
       IF(EQREAL(PNUMS(17),0.))THEN                                  
         IREFR = 0                                                     
         CALL MSGERR(0, 'Refraction deactivated')                      
       ENDIF
     ENDIF
!
!    *** numerical scheme in frequency space :                  ***
!    *** 1) Fully implicit scheme in frequency space (iterative ***
!    ***    SIP solver)                                         ***
!    *** 2) Explicit scheme in frequency space:                 ***
!    ***    Energy is removed from the spectrum based on a CFL  ***
!    ***    criterion                                           ***
!    *** 3) Explicit scheme in frequency space:                 ***
!    ***    No CFL limitation near blocking point --> unstable  ***
!    ***    integration. Therefore a filter is applied with a   ***
!    ***    diffusion coeff.                                    ***
!
     CALL INCSTR('SIGIMPL',SIGIM,'UDF',' ')

     IF(SIGIM == 'SIGIM')THEN                     
!      *** implicit solver ***
!      ***  This is the default option   PNUMS(8) = 1.  ***
       PNUMS(8) = 1.
       CALL INREAL ('CSS'   , PNUMS(7), 'UNC', 0.)
       CALL INREAL ('EPS1'  , PNUMS(11) , 'UNC', 0.)
       CALL INREAL ('EPS2'  , PNUMS(12) , 'UNC', 0.)
       CALL INREAL ('OUTP'  , PNUMS(13) , 'UNC', 0.)
       CALL INREAL ('NITER' , PNUMS(14), 'UNC', 0.)
     ELSE IF(SIGIM == 'SIGEX')THEN                 
!
!      *** 2) explicit scheme ***
       PNUMS(8) = 2.
       CALL INREAL ('CFL'  , PNUMS(19), 'UNC', 0.)
!
     ELSE IF(SIGIM == 'FIL')THEN
!      *** 3) explicit scheme -> filter the spectrum  ***
       PNUMS(8) = 3.
       CALL INREAL ('DIFFC'  , PNUMS(9), 'UNC', 0.)
!
     END IF

     SETUP = .FALSE.
     CALL INLOGC('SETUP',SETUP,'DFL','T')

     IF(SETUP)THEN                                         
!      *** iterative solver        ***                                 
!      *** Settings for the setup  ***                                 
       CALL INREAL ('EPS2' , PNUMS(23), 'UNC', 0.)                     
       CALL INREAL ('OUTP' , PNUMS(24), 'UNC', 0.)                     
       CALL INREAL ('NITER', PNUMS(25), 'UNC', 0.)                     
     ENDIF                                                             
   ENDIF
!
!     ------------------------------------------------------------------
!
!     **  Command COMPUTE   **                                       
!
!       ==============================================================
!
!                  |  STATionary  [time]                      |
!       COMPute ( <                                            > )        
!                  |                    | -> Sec  |           |
!                  |  ([tbegc] [deltc] <     MIn   > [tendc]) |
!                                       |    HR   |
!                                       |    DAy  |
!
!       ==============================================================
!

   CALL INCSTR('COMPUT',COMPUT,'REQ',' ')
   IF(COMPUT == 'COMP')THEN
     RUNMADE = .TRUE.                                                  
!     CALL INLOGC('COM_STAT',COM_STAT,'REQ','T')
     COM_STAT = .FALSE.
     IF(NSTATM <= 0 .OR. COM_STAT)THEN                         
       IF(NSTATM == -1) NSTATM = 0                                    
       IF(NSTATM > 0) CALL INCTIM (ITMOPT,'TIME',TINIC,'REQ',0.)     
       IF(TINIC < TIMCO)THEN                                      
         CALL MSGERR (2, '[time] before current time')                 
         TINIC = TIMCO                                                 
       ENDIF                                                           
       TFINC = TINIC                                                   
       TIMCO = TINIC                                                   
       DTW = 1.E10                                                      
       RDTIM = 0.                                                      
       NSTATC = 0                                                      
       MTC = 1                                                         
     ELSE

       TINIC = IFLBEG(IGRD)
       IF(TINIC < TIMCO)THEN                                      
         CALL MSGERR (2, 'start time [tbegc] before current time')     
         TINIC = TIMCO                                                 
       ENDIF                                                           
       CALL INREAL('NS_DELTC',DTW,'REQ',1)
       
       CALL INCSTR('TIME_UNIT',UNIT,'UNC','SECOND')
       IF(TRIM(UNIT) == 'DAY' .OR. TRIM(UNIT) == 'day')THEN
         DTW = DTW*24.0*3600.0
       ELSE IF(TRIM(UNIT) == 'HOUR' .OR. TRIM(UNIT) == 'hour')THEN
         DTW = DTW*3600.0
       ELSE IF(TRIM(UNIT) == 'MINUTE' .OR. TRIM(UNIT) == 'minute')THEN
         DTW = DTW*60.0
       ELSE IF(TRIM(UNIT) /= 'SECOND' .AND. TRIM(UNIT) /= 'second')THEN
         CALL MSGERR(3, 'the unit of time step must be DAY,HOUR,MINUTE &
	                 or SECOND')	 
       END IF	 	 	 
       
       NSTATC = 1 
       IMDTW = SECONDS2TIME(DTW)
       TFINC = SECONDS(EndTime)                                                     
!
!      *** tfinc must be greater than tinic **
       DIFF = TFINC - TINIC
       
       if(msr) PRINT*,'TINIC=',TINIC,'TFINC=',TFINC,'DIFF=',DIFF,'MYID=',MYID
       IF(DIFF <= 0.) CALL MSGERR (3,                                &
              'start time [tbegc] greater or equal end time [tendc]')
!
!        **The number of computational steps is calculated
       RDTIM = 1./DTW
       MTC = NINT ((TFINC - TINIC)/DTW)                                 
       IF(MOD(TINIC-TFINC,DTW) > 0.01*DTW .AND.                        &
          MOD(TINIC-TFINC,DTW) < 0.99*DTW)                             &
	       CALL MSGERR (1,                                       &
	       'DTW is not a fraction of the computational period')
       TIMCO = TINIC
     ENDIF
!JQI     IF(NSTATM > 0) CHTIME = DTTIWR(ITMOPT, TIMCO)                   
     NCOMPT = NCOMPT + 1                                               
!JQI     IF(NCOMPT > 50) CALL MSGERR (2,                                 &
!JQI                   'No more than 50 COMPUTE commands are allowed')  
     RCOMPT(NCOMPT,1) = REAL(NSTATC)                                   
     RCOMPT(NCOMPT,2) = REAL(MTC)                                      
     RCOMPT(NCOMPT,3) = TFINC                                          
     RCOMPT(NCOMPT,4) = TINIC                                          
     RCOMPT(NCOMPT,5) = DTW                                             
!    set ITERMX equal to MXITST in case of stationary computations
!    and to MXITNS otherwise
     IF(NSTATC == 0)THEN
       ITERMX = MXITST                                                 
     ELSE
       ITERMX = MXITNS                                                 
     ENDIF
!     RETURN                                                            
   ENDIF

   IF(MSR)PRINT*,'AT THE END OF SWREAD, MCGRD = ',MCGRD
!   STOP
   RETURN
   END SUBROUTINE SWREAD
 
!***********************************************************************
!                                                                      *
   SUBROUTINE SSFILL (SPCSIG, SPCDIR)                                  
!                                                                      *
!***********************************************************************
!
!     Discretisation in frequency (sigma) and direction (theta)
!
!  Method
!
!     Create a logarithmic distribution in frequency between lowest and highest
!     frequencies and store them in the array SPCSIG.
!
!     Create a distribution in direction that does not coincide with the orientation
!     of the computational grid, calculate some geometric derivatives and store
!     it in the array SPCDIR
!
!***********************************************************************
!
   USE OCPCOMM1                                                        
   USE OCPCOMM2                                                        
   USE OCPCOMM3                                                        
   USE OCPCOMM4                                                        
   USE SWCOMM1                                                         
   USE SWCOMM2                                                         
   USE SWCOMM3                                                         
   USE SWCOMM4                                                         
   USE M_WCAP  
      
   IMPLICIT NONE                                                        

   REAL    SPCDIR(MDC,6)                                               
   REAL    SPCSIG(MSC)     
   
   REAL    :: SFAC,OLDDIR
   INTEGER :: ID,IS                                            

!
!  Allocate arrays in module M_WCAP                                    
!
   IF (.NOT.ALLOCATED(SIGPOW)) ALLOCATE (SIGPOW(MSC,6))                
!
!
!  distribution of spectral frequencies
!
!  FRINTF is the frequency integration factor (=df/f)                  
   FRINTF = ALOG(SHIG/SLOW) / FLOAT(MSC-1)                             
   SFAC   = EXP(FRINTF)                                                
   FRINTH = SQRT(SFAC)
!  determine spectral frequencies (logarithmic distribution)
   SPCSIG(1) = SLOW                                                    
   DO IS = 2, MSC
     SPCSIG(IS) = SPCSIG(IS-1) * SFAC                                 
   END DO  
!
!  Calculate powers of sigma and store in global array
!
   SIGPOW(:,1) = SPCSIG                                                
   SIGPOW(:,2) = SPCSIG**2                                             
   SIGPOW(:,3) = SPCSIG * SIGPOW(:,2)                                  
   SIGPOW(:,4) = SPCSIG * SIGPOW(:,3)                                  
   SIGPOW(:,5) = SPCSIG * SIGPOW(:,4)                                  
   SIGPOW(:,6) = SPCSIG * SIGPOW(:,5)                                  
!
!  distribution of spectral directions
!
   DO ID = 1, MDC
     SPCDIR(ID,1) = SPDIR1 + FLOAT(ID-1)*DDIR                         
!    if a direction coincides with a direction of the (regular)       
!    computational grid it is slightly changed
     IF(OPTG == 1)THEN                                              
       IF(ABS(MODULO(SPCDIR(ID,1)-ALPC+0.25*PI_W ,0.5*PI_W)-0.25*PI_W)  &
	      < 1.E-6)THEN
         OLDDIR = SPCDIR(ID,1)
         SPCDIR(ID,1) = OLDDIR + 2.E-6                                
       ENDIF
     ENDIF                                                            
     SPCDIR(ID,2) = COS(SPCDIR(ID,1))                                 
     SPCDIR(ID,3) = SIN(SPCDIR(ID,1))                                 
     SPCDIR(ID,4) = SPCDIR(ID,2) **2                                  
     SPCDIR(ID,5) = SPCDIR(ID,2) * SPCDIR(ID,3)                       
     SPCDIR(ID,6) = SPCDIR(ID,3) **2                                  
   END DO  
!
   RETURN
   END SUBROUTINE SSFILL
