SUBROUTINE wamodel (nadv, nsubadv)
! ----------------------------------------------------------------------
!**** *WAMODEL* - 3-G WAM MODEL - TIME INTEGRATION OF WAVE FIELDS.
!     S.D. HASSELMANN  MPI       1.12.85
!     G. KOMEN         KNMI         6.86  MODIFIED FOR SHALLOW WATER
!     P. JANSSEN                          ASPECTS.
!     S.D. HASSELMANN  MPI       15.2.87  MODIFIED FOR CYBER 205.
!     P. LIONELLO      ISDGM      6.3.87  MODIFIED TO OUTPUT SWELL.
!     S.D. HASSELMANN  MPI        1.6.87  ALL VERSIONS COMBINED INTO
!                                         ONE MODEL. DEEP AND SHALLOW
!                                         WATER , CRAY AND CYBER 205
!                                         VERSION.
!     CYCLE_2 MODICIFATIONS:
!     ----------------------
!     L. ZAMBRESKY     GKSS        10.87  OPTIMIZED FOR CRAY, CYBER 205
!     H. GUNTHER
!     A. SPEIDEL       MPI          4.88  VARIABLE DIMENSIONS, INTERNAL
!                                         CHECKS (CFL-CRITERION).
!     A. SPEIDEL       MPI         11.88  CHANGES FOR CRAY-2.
!     K. HUBBERT       POL          6.89  DEPTH AND CURRENT REFRACTION.
!                                         PRECALCULATION OF TERMS IN
!                                         *PROPDOT*.
!                                         SOLVE WAVE ACTION EQUATION
!                                         FOR CURRENT REFRACTION.
!     CYCLE_3 MODICIFATIONS:
!     ----------------------
!     R. PORTZ , S.D. HASSELMANN   MPI          1990
!      - RESTRUCTURE MODEL TO CALL THE ACTUAL INTEGRATION IN TIME
!        AS A SUBROUTINE: WAMODEL. A SHELL PROGRAM "WAMSHELL" READS
!        OUTPUT FROM PREPROC AND COMPUTES THE WIND ARRAYS FOR THE
!        INTEGRATION PERIOD FROM PREWIND, WHICH HAS BEEN INCORPORATED
!        AS A SUBROUTINE.
!      - ALL INTERMEDIATE AND RESTART I/O IS DONE IN THE SUBROUTINE
!        WAMODEL AND INPREST.
!      - THE COMMON BLOCK IN THE PREPROCESSOR AND MODEL ARE MADE
!        COMPATIBLE.
!      - THE COMPUTATION OF SEVERAL PARAMETERS HAS BEEN TRANSFERRED
!        FROM THE MODEL TO PREPROC.
!      - DEPTH AND CURRENT REFRACTION HAS BEEN INCORPORATED INTO THE
!        MODEL.
!      - OPEN BOUNDARIES ARE INCORPORATED IN THE MODEL.
!      - SEVERAL MINOR ERRORS HAVE BEEN REMOVED.
!      - THE BUFFERED I/O FOR THE CYBER 205 HAS BEEN CHANGED INTO A
!        BINARY READ AND WRITE.
!     CYCLE_4 MODICIFATIONS:
!     ----------------------
!     L. ZAMBRESKY   GKSS/ECMWF   6.89  ECMWF SUB VERSION
!                                       BASED ON CYCLE_2.
!     H. GUNTHER     GKSS/ECMWF 10.89  ECMWF SUB VERSION REORGANIZED.
!                                      - COMMON BLOCK STRUCTURE.
!                                      - BLOCKING STRUCTURE.
!                                      - TIME COUNTING.
!                                      - GRIDDED OUTPUT FIELDS.
!                                      - HEADERS ADDED TO OUTPUT FILES.
!                                      - ERRORS IN PROPAGATION CORRECTED
!     P.A.E.M. JANSSEN KNMI      1990  COUPLED MODEL.
!     H. GUNTHER     GKSS/ECMWF  8.91  LOGARITHMIC DEPTH TABLES.
!                                      MPI CYCLE_3 AND ECMWF VERSIONS
!                                      COMBINED INTO CYCLE_4.
!SHALLOW
!          DIFFERENCES FOR SHALLOW WATER RUNS TO DEEP WATER RUNS
!          ARE ENCLOSED IN COMMENT LINES : 'CSHALLOW'.
!SHALLOW
!NEST
!          DIFFERENCES FOR NESTED GRID RUNS TO NORMAL RUNS
!          ARE ENCLOSED IN COMMENT LINES : 'CNEST'.
!NEST
!REFRA
!          DIFFERENCES FOR REFRACTION RUNS TO NORMAL RUNS
!          ARE ENCLOSED IN COMMENT LINES : 'CREFRA'.
!REFRA
!*    PURPOSE.
!     --------
!       COMPUTATION OF THE 2-D FREQUENCY-DIRECTION WAVE SPECTRUM AT ALL
!       GRID POINTS FOR A GIVEN INITIAL SPECTRUM AND FORCING SURFACE
!       STRESS FIELD.
!**   INTERFACE.
!     ----------
!       *CALL* *WAMODEL (NADV)*
!         *NADV*     INTEGER   NUMBER OF ADVECTION ITERATIONS.
!     METHOD.
!     -------
!       GRID POINTS ARE LAT - LONG,VECTORIZATION IS ACHIEVED BY RUNNING
!       THROUGH THE GRID POINTS IN AN INNER LOOP ORGANIZED AS 1-D ARRAY
!       IN BLOCKS,-ALL COMPUTATIONS ARE CARRIED OUT FOR ONE BLOCK AT A
!       TIME (SEE "BLOCK STRUCTURE" BELOW)
!       ALL COMPONENTS OF THE SPECTRUM ARE COMPUTED PROGNOSTICALLY FROM
!       THE SPECTRAL TRANSPORT EQUATION UP TO A VARIABLE CUT-OFF
!       FREQUENCY = MAX(4*FPM,2.5*FMEAN),WHERE FPM IS THE
!       PIERSON MOSKOVITZ FREQUENCY AND FMEAN IS THE MEAN FREQUENCY,
!       BEYOND THE PROGNOSTIC CUTOFF A DIAGNOSTIC F**-5 TAIL IS ATTACHED
!       CONTINUOUSLY FOR EACH DIRECTION,
!       SOURCE FUNCTIONS ARE TAKEN FROM KOMEN ET AL(1984)
!       THE NONLINEAR TRANSFER IS PARAMETERIZED BY THE DISCRETE INTER-
!       ACTION APPROXIMATION OF HASSELMANN ET AL (1985B)
!       THE SOURCE FUNCTION AND THE ADVECTION TERM ARE INTEGRATED ON TWO
!       DIFFERENT TIME STEP LEVELS AND WITH DIFFERENT METHODS,-THE
!       ADVECTION TIME STEP IS A MULTIPLE OF THE SOURCE FUNCTION
!       TIME STEP.
!       THE SOURCE FUNCTIONS ARE INTEGRATED IMPLICITLY ACCORDING TO
!       HASSELMANN AND HASSELMANN (1985A),-THE RELEVANT FUNCTIONAL
!       DERIVATIVES OF THE INDIVIDUAL SOURCE FUNCTIONS REQUIRED FOR THE
!       SOLUTION OF THE IMPLICIT EQUATION ARE COMPUTED WITHIN THE SOURCE
!       FUNCTION SUBS,- THE TIME STEP IS TYPICALLY 20 MIN,
!       THE ADVECTION IS INTEGRATED BY A FIRST ORDER UPWIND SCHEME,ALSO
!       ACCORDING TO HASSELMANN AND HASSELMANN (1985A),-THE ADVECTIVE
!       TIMESTEP IS DEPENDENT ON THE FREQUENCY AND SPATIAL GRID IN
!       ACCORDANCE WITH CFL,
!       WINDS ARE READ IN EVERY WIND TIME STEP.IF THE WIND TIME STEP IS
!       GREATER THAN THE SOURCE TERM TIME STEP DELTWIND/DELTSOURCE STEPS
!       ARE INTEGRATED WITH CONSTANT WINDS,
!       WIND TIME STEP,PROPAGATION TIME STEP AND SOURCE TERM TIME STEP
!       SHOULD HAVE INTEGER RATIOS, THEY ARE GIVEN IN SECONDS AT
!       FULL MINUTES.
!NEST
!       ZERO ENERGY INFLUX IS ASSUMED AT COAST LINES. OPEN BOUNDARIES
!       ARE INCORPORATED IN THE MODEL, IF IT RUNS AS A NESTED GRID.
!NEST
!       BLOCK STRUCTURE (SEE PREPROC FOR DETAILS):
!       SEA POINTS ARE COLLECTED INTO A 1-DIMENSIONAL ARRAY.
!       BLOCKS OF MAXIMALLY NIBLO ELEMENTS.
!       SEA POINTS ARE COUNTED ALONG LINES OF LATITUDES FROM LEFT COAST
!       TO RIGHT COAST WORKING FROM SOUTH TO NORTH.
!       BLOCKS OVERLAP OVER TWO LATITUDE LINES,TO COMPUTE NORTH-SOUTH
!       ADVECTION TERMS, SEE ALSO COMMON GRIDPAR AND UBUF.
!       THE WIND FILES FOR THE BLOCKED WINDS CREATED BY PREWIND ARE
!       READ AND DELETED IN SUB IMPLSCH (IU17 AND IU18). THE FILE
!       NAMES ARE CREATED IN SUB CREWFN AND AN IMPLICIT OPEN IS USED.
!       THE FILE HANDLING SUBS OPENFIL, GSFILE AND CREWFN ARE COMPUTER
!       DEPENDENT AND MAY BE ADOPTED BY THE USER.
!       THE PROGRAM CLOSES AND DELETES ALL WORK FILES.
!       ALL PARAMETERS HAVE TO BE THE VALUES GIVEN AT THE END OF THE
!       PREPROC OUTPUT IN COLUMN 'REQUIRED'.
!     EXTERNALS.
!     ----------
!       *ABORT*     - TERMINATES PROCESSING.
!       *AIRSEA*    - SURFACE LAYER STRESS.
!NEST
!       *BOUINPT*   - BOUNDARY VALUE INPUT.
!NEST
!       *CREWFN*    - CREATES A WIND FILE NAME.
!REFRA
!       *DOTDC*     - READ COMMON REFDOT.
!REFRA
!       *FEMEAN*    - COMPUTATION OF MEAN FREQUENCY AT EACH GRID POINT.
!       *FILLBL*    - ADD LATITUDES TO A BLOCK.
!       *GSFILE*    - ROUTINE TO DYNAMICALLY FETCH OR DISPOSE FILES.
!       *HEADBC*    - WRITE BOUNDARY OUTPUT FILE HEADER.
!       *iecf_len*  - LENGTH OF A CHARACTER ARRAY.
!       *IMPLSCH*   - IMPLICIT SCHEME FOR INTEGRATION OF SOURCE
!                     FUNCTIONS IN TIME AND INPUT OF WINDS.
!       *INCDATE*   - UPDATE DATE TIME GROUP.
!REFRA
!       *INTPOL*    - MAP SPECTRUM FROM SIGMA TO OMEGA SPACE.
!REFRA
!NEST
!       *INTSPEC*   - INTERPOLATION OF SPECTRA.
!NEST
!       *MAKEGRID*  - MAKES GRIDDED FIELDS.
!       *OPENFIL*   - OPENS A FILE.
!       *OUTBC*     - OUTPUT OF BOUNDARY VALUES.
!       *OUTBS*     - CONTROLS OUTPUT FROM BLOCKS.
!       *OUTGRID*   - SAVE BLOCKED PARAMETERS INTO GRID ARRAYS.
!       *OUTINT*    - OUTPUT OF INTEGRATED PARAMETERS.
!       *OUTPP*     - WRITES GRIDDED FIELDS TO PRINTER.
!       *OUTSPP*    - OUTPUT OF SPECTRA AT SELECTED POINTS.
!       *PEAKFR*    - COMUTE PEAK FREQUENCY.
!       *PROPAG*    - PROPAGATION SCHEME.
!       *PRSPP*     - PRINT A SPECTRUM.
!       *PRSPPS*    - PRINT A SWELL SPECTRUM.
!NEST
!       *ROTSPEC*   - ROTATE A SPECTRUM.
!NEST
!       *SAVREST*   - DISPOSE RESTART FILES.
!SHALLOW
!       *SBOTTOM*   - COMPUTES BOTTOM DISSIPATION SOURCE TERM AND
!                     LINEAR CONTRIBUTION TO FUNCTIONAL MATRIX.
!SHALLOW
!       *SDISSIP*   - COMPUTATION OF DISSIPATION SOURCE FUNCTION
!                     AND LINEAR CONTRIBUTION OF DISSIPATION TO
!                     FUNCTIONAL MATRIX IN IMPLICIT SCHEME.
!       *SEMEAN*    - COMPUTATION OF TOTAL ENERGY AT EACH GRID POINT.
!       *SEPWISW*   - COMPUTATION OF 2-DIMENSIONAL SWELL DISTRIBUTION
!                     TOTAL SWELL ENERGY, MEAN SWELL DIRECTION, AND
!                     MEAN SWELL FREQUENCY AT EACH GRID POINT.
!       *SINPUT*    - COMPUTATION OF INPUT SOURCE FUNCTION, AND
!                     LINEAR CONTRIBUTION OF INPUT SOURCE FUNCTION
!                     TO FUNCTIONAL MATRIX IN IMPLICIT SCHEME.
!       *SNONLIN*   - COMPUTATION OF NONLINEAR TRANSFER RATE AND
!                     DIAGONAL LINEAR CONTRIBUTION OF NONLINEAR SOURCE
!                     FUNCTION TO FUNCTIONAL MATRIX.
!       *SPLITBL*   - STORING OVERLAPPING LATITUDES OF BLOCKS.
!       *STHQ*      - COMPUTATION OF MEAN WAVE DIRECTION AT EACH
!                     GRID POINT.
!       *STRESSO*   - COMPUTATION OF WAVE STRESS.
!NEST
!       *STRSPEC*   - STRETCH A SPECTRUM.
!NEST
!     REFERENCE.
!     ----------
!       SNYDER, R.L., F.W. DOBSON, J.A. ELLIOT, AND R.B. LONG:
!          ARRAY MEASUREMENTS OF ATMOSPHERIC PRESSURE FLUCTUATIONS
!          ABOVE SURFACE GRAVITY WAVES. J.FLUID MECH. 102, 1-59 ,1981.
!       G. KOMEN, S. HASSELMANN, K. HASSELMANN:
!          ON THE EXISTENCE OF A FULLY DEVELOPED WIND SEA SPECTRUM.
!          JPO,1984.
!       S. HASSELMANN, K. HASSELMANN, J.H. ALLENDER, T.P. BARNETT:
!          IMPROVED METHODS OF COMPUTING AND PARAMETERIZING THE
!          NONLINEAR ENERGY TRANSFER IN A GRAVITY WAVE SPECTRUM.
!          JPO, 1985.
!       S. HASSELMANN, K. HASSELMANN: A GLOBAL WAVE MODEL,
!          WAM REPORT,JUNE,30/1985.
!       P. JANSSEN, G. KOMEN: A SHALLOW WATER EXTENSION OF THE
!          3-G WAM-MODEL. WAM REPORT 1985.
!       THE WAMDI GROUP: THE WAM MODEL - A THIRD GENERATION OCEAN
!          WAVE PREDICTION MODEL. JPO, VOL. 18, NO. 12, 1988.
!       P.A.E.M JANSSEN: JPO, 1989 AND 1991.
!       K. HASSELMANN: TRANSPORT EQUATION OF FINITE DEPTH SURFACE
!          WAVE SPECTRUM IN TIME DPENDANT CURRENT AND DEPTH FIELD USING
!          NONCANONICAL SPACIAL (SPHERICAL) AND WAVE NUMBER (FRQUENCY-
!          DIRECTION) COORDINATES. WAM REPROT 1988.
! ----------------------------------------------------------------------
USE params  ; USE coupl  ; USE fredir ; USE gridpar ; USE indnl   ; USE intpar
USE intspar ; USE map    ; USE meanpa ; USE ndtf    ; USE over    ; USE cout
USE shallow ; USE source ; USE spe    ; USE idfout  ; USE idfout1 ; USE stat
USE swell   ; USE table  ; USE testo  ; USE text    ; USE ubuf    ; USE units
USE wind    ; USE couple_mpi
!NEST
USE cbound ; USE fbound
!NEST
!REFRA
USE current ; USE refdot
!REFRA
#ifdef ice
USE ice1
#endif
#ifdef linux
USE flushfiles
#endif
#ifndef _OPENMP
USE my_omp
#endif
IMPLICIT NONE

INTEGER, INTENT(IN) :: nadv, nsubadv
#ifdef ice
INTEGER   :: mask_ice(ngx,ngy)
#endif
CHARACTER (LEN=12) :: idtpra, idrn, idtprf
INTEGER   :: out_interval
CHARACTER :: srcmd * 88, rest * 11, mk_dir * 51
INTEGER   :: j, k, m, id, ig, ij, iubc, ifil, kadv, nest
INTEGER   :: it1, it2, iareas, isub, ijss, ijll, ijsst, ijllt, igh1
REAL      :: xd
#ifdef _OPENMP
INTEGER   :: omp_get_thread_num
#endif

! Interval (in sec) for restart files saving
out_interval = 43200

#ifdef ice
OPEN(66, FILE = 'ICEMASK_FILE', FORM = 'UNFORMATTED')
READ(66) mask_ice
CLOSE(66)
WRITE(iu06,*) 'ICE MASK'
DO i = 1, ngx
   WRITE(iu06, '(157i1)') (mask_ice(i,j), j = 1, ngy)
ENDDO
DO ig = 1, igl
   CALL trans_grid (icemask, mask_ice, ijs(ig), ijl(ig), ig)
ENDDO
WRITE(iu06,*) '   ICE MASK READ AND TRANSFORMED '
#endif
!*    1. ADVECTION TIME LOOP.
!        --------------------
DO kadv = 1, nadv
!*    1.1 FIX END DATE OF THIS PROPAGATION STEP AND OUTPUT TIMES.
!         -------------------------------------------------------
   idtpra = idtpro
   CALL incdate (idtpro, idelpro)
   IF( itest .GE. 2 ) THEN
      WRITE(iu06,*) '   SUB. WAMODEL: START OF PROPAGATION '
      WRITE(iu06,*) '     START DATE IS    IDTPRA = ', IDTPRA
      WRITE(iu06,*) '     END DATE WILL BE IDTPRO = ', IDTPRO
   ENDIF
!         UPDATE OUTPUT TIMES.
   IF( noutt .GT. 0 ) THEN
      idtintt = '000000000000'
      idtints = '000000000000'
      idtspt  = '000000000000'
      idtsps  = '000000000000'
      DO j = 1, noutt
         IF( idtpro .EQ. ioutt(j) ) THEN
            IF( pflag20 .OR. fflag20 )   idtintt = ioutt(j)
            IF( pflag21 .OR. fflag21 )   idtints = ioutt(j)
            IF( pflag25 .OR. fflag25 )   idtspt  = ioutt(j)
            IF( pflag26 .OR. fflag26 )   idtsps  = ioutt(j)
         ENDIF
      ENDDO
   ELSE
      IF((pflag20 .OR. fflag20) .AND. idtintt .LT. idtpro)   CALL incdate (idtintt,idelint)
      IF((pflag21 .OR. fflag21) .AND. idtints .LT. idtpro)   CALL incdate (idtints,idelins)
      IF((pflag25 .OR. fflag25) .AND. idtspt  .LT. idtpro)   CALL incdate (idtspt,idelspt)
      IF((pflag26 .OR. fflag26) .AND. idtsps  .LT. idtpro)   CALL incdate (idtsps,idelsps)
   ENDIF
   IF( idtres .LT. idtpro )     CALL incdate (idtres,idelres)
   IF( idtout .LT. idtpro )     CALL incdate (idtout,out_interval)
!*    1.2 CHECK FILE STATUS OF BUFFER OPERATIONS AND RESET FILES.
!         ------------------------------------------------------_
   REWIND iu12
   IF( igl .NE. 1 ) THEN
      REWIND iu08
      REWIND iu11
      REWIND iu13
      REWIND iu14
!*    1.3 CHANGE WORKING FILE UNITS.
!         --------------------------
      ifil = iu11
      iu11 = iu12
      iu12 = ifil
      ifil = iu13
      iu13 = iu14
      iu14 = ifil
!*    1.4 READ FIRST BLOCK OF SPECTRA, IF MULTI BLOCK VERSION.
!         ----------------------------------------------------
      READ (iu11) (((fl1(ij, k, m), ij = 0, niblo), k = 1, nang), m = 1, nfre)
   ENDIF
!*    1.5  LOOP FOR BLOCKS OF LATITUDES.
!          -----------------------------
   DO ig = 1, igl
!*    1.5.1 EVEN OR ODD BLOCK NUMBER
!          -------------------------
      igh1 = (ig / 2) * 2
!*    1.5.2 INPUT NEIGHBOURING GRID POINT INDICES (COMMON BLOCK UBUF).
!           ----------------------------------------------------------
      IF( igl .NE. 1 ) THEN
         READ (iu08) klat, klon
!SHALLOW
!*    1.5.3 COMPUTE SHALLOW WATER TABLE INDICES.
!           ------------------------------------
         IF( ishallo .NE. 1 ) THEN
            DO ij = 1, ijlt(ig)
               xd        = LOG (depth(ij,ig) / deptha) / LOG (depthd) + 1.
               id        = NINT (xd)
               id        = MAX (id, 1)
               indep(ij) = MIN (id, ndepth)
            ENDDO
         ENDIF
!SHALLOW
!*    1.5.4 COUPLING WITH NEIGHBOURING BLOCKS IG +- 1 AND START
!*          INPUT OF SPECTRA FOR BLOCK IG+1.
!           ----------------------------------------------------
         IF( igh1 .EQ. ig ) THEN
            CALL fillbl (fl2, ig, iu13)
            IF( ig .NE. igl )   READ (iu11) (((fl1(ij,k,m),ij=0,niblo),k=1,nang),m=1,nfre)
         ELSE
            CALL fillbl (fl1, ig, iu13)
            IF( ig .NE. igl )   READ (iu11) (((fl2(ij,k,m),ij=0,niblo),k=1,nang),m=1,nfre)
         ENDIF
      ENDIF

!$OMP PARALLEL private(isub)
#ifdef ice
!        SET FL1 ON ICE POINTS TO ZERO
!NNIKOS I thing this is how it should be for the parallel version to work correctly in any case (even if ig>1)
      IF(ig .EQ. 1) THEN
         CALL setice (fl1, ITASK(1, OMP_GET_THREAD_NUM()+1),     &
                           ITASK(2, OMP_GET_THREAD_NUM()+1), ig)
      ELSE
         CALL setice (fl1, ijs(ig), ijl(ig), ig)
      ENDIF
#endif
!*    1.5.5 COMPUTATION OF PROPAGATION.
!           ---------------------------
!           PROPAGS1     QUADRANT CO-ORDINATE
!           PROPAGS2     OCTANT CO-ORDINATE (ONLY WITH IREFRA = 0, 1)
      DO isub = 1, nsubadv
#ifdef propagation2
!NNIKOS should be ijsst/ijllt (or itask(3,)/itask(4,) instead of 1,niblo
         CALL propags2 (isub, ITASK(1, OMP_GET_THREAD_NUM()+1),          &
                              ITASK(2, OMP_GET_THREAD_NUM()+1), 1, niblo)
#elif propagation1
!NNIKOS change propags1 as well.
         CALL propags1 (isub, ITASK(1, OMP_GET_THREAD_NUM()+1),          &
                              ITASK(2, OMP_GET_THREAD_NUM()+1), 1, niblo)
#endif
!wluipro-end
      ENDDO
!$OMP END PARALLEL
!*    1.5.6 INTEGRATION OF SOURCE TERMS OVER SUB TIME STEPS BETWEEN
!*          PROPAGATION TIME STEPS.
!           -------------------------------------------------------
      CALL implsch (fl3, fl1, diss,ijs(ig), ijl(ig), ig, igl, idtpra)
!NEST
#ifdef ice
!  SET FL3 ON ICE POINTS TO ZERO
!NNIKOS I thing this is how it should be for the parallel version to work correctly in any case (even if ig>1)
!$OMP PARALLEL
      IF(ig .EQ. 1) THEN
         CALL setice (fl3, ITASK(1, OMP_GET_THREAD_NUM()+1),     &
                           ITASK(2, OMP_GET_THREAD_NUM()+1), ig)
      ELSE
         CALL setice (fl3, ijs(ig), ijl(ig), ig)
      ENDIF
!$OMP END PARALLEL
#endif
!*    1.5.7 INPUT OF BOUNDARY VALUES.
!           -------------------------
      IF( ibounf .EQ. 1 ) THEN
         CALL bouinpt (fl3, ig, iu02)
      ENDIF
!NEST
!*    1.5.8 OUTPUT OF SECOND LATITUDES FOR NEXT TIMESTEP
!*          START OF OUTPUT OF ALL SPECTRA IN BLOCK.
!           --------------------------------------------
      IF( igl .NE. 1 ) THEN
         IF( ig .NE. 1 )   CALL splitbl (fl3, ijs(ig), ijl2(ig), iu14)
         WRITE(iu12) (((fl3(ij,k,m), ij=0, niblo), k=1, nang), m=1, nfre)
      ENDIF
!NEST
!*    1.5.9 OUTPUT OF BOUNDARY POINTS.
!           --------------------------
      IF( ibounc .EQ. 1 ) THEN
         iubc = 50
         DO nest = 1, nests
            CALL outbc (fl3, ijs(ig), ijl(ig), ig, iubc, nest)
            iubc = iubc + 1
         ENDDO
      ENDIF
!NEST
!*    1.5.10 MODEL OUTPUT INTEGRATED DATA ARE SAVED IN COMMON.
!*           SPECTRA ARE SAVED ON DISK AND PRINTED.
!            -------------------------------------------------
!NIKOS Spaces are placed according to MY logical perception
      IF(     idtintt .EQ. idtpro  .OR.  idtints .EQ. idtpro  &
         .OR. idtspt  .EQ. idtpro  .OR.  idtsps  .EQ. idtpro)   THEN
         IF( igh1 .EQ. ig ) THEN
!REFRA
            IF( irefra .EQ. 2 ) THEN
               CALL intpol (fl3, sl, ijs(ig), ijl(ig), ig)
               CALL outbs (sl, fl2, ijs(ig), ijl(ig), ig, iu25, iu26)
!REFRA
            ELSE
               CALL outbs (fl3, fl2, ijs(ig), ijl(ig), ig, iu25, iu26)
            ENDIF
         ELSE
!REFRA
            IF( irefra .EQ. 2 ) THEN
               CALL intpol (fl3, sl, ijs(ig), ijl(ig), ig)
               CALL outbs (sl, fl1, ijs(ig), ijl(ig), ig, iu25, iu26)
!REFRA
            ELSE
               CALL outbs (fl3, fl1, ijs(ig), ijl(ig), ig, iu25, iu26)
            ENDIF
         ENDIF
      ENDIF
!*    BRANCHING BACK TO 1.5 FOR NEXT BLOCK OF LATITUDES
   ENDDO
!*    1.6 IF ONE BLOCK VERSION COPY RESULTS.
!         ----------------------------------
   IF( igl .EQ. 1 ) THEN
      DO m = 1, nfre
         DO k = 1, nang
!$OMP PARALLEL DO schedule(static)
            DO ij = 0, niblo
               fl1(ij, k, m) = fl3(ij, k, m)
            ENDDO
!$OMP END PARALLEL DO
         ENDDO
      ENDDO
   ENDIF
!*    1.7 ONE PROPAGATION TIMESTEP DONE FOR ALL BLOCKS.
!         ---------------------------------------------
!         WRITE INTEGRATED DATA TO FILE AND/OR PRINTER
!         DATA WERE COLLECTED INSIDE THE BLOCK LOOP.
   IF( idtintt.EQ.idtpro .OR. idtints.EQ.idtpro )   CALL outint (iu06, iu20, iu21)
!         PRINT TIME.
   WRITE(iu06, 112) idtpro
  112 FORMAT(/,3X,'!!!!!!!!!!!!!! ', 'WAVE FIELDS INTEGRATED  DATE IS: ',&
             a10, '  !!!!!!!!!!!!!! ')
!*    1.8 OUTPUT FILES AND RECOVERY FILES ARE DISPOSED WHEN
!         TIME REACHES THE DISPOSE DATE.
!         -------------------------------------------------
   IF( idtres.EQ.idtpro .OR. idatee.EQ.idtpro ) THEN
!*    1.8.1 DATE FOR NEW RESULT FILES TO BE OPENDED.
!           ----------------------------------------
      IF( idtpro .LT. idatee ) THEN
         idrn = idtres
         CALL incdate (idrn, idelres)
      ELSE
         idrn = '000000000000'
      ENDIF
!*    1.8.2 STORE 2-D SPECTRA AT SELECTED GRID POINTS,
!           INTEGRATED PARAMETERS OF ENTIRE GRID.
!           ------------------------------------------
!           SAVE INTEGRATED DATA FILE.
      IF(fflag20)     CALL gsfile (iu06, iu20, 1, idrn, 'MAP', 'S')
!           SAVE INTEGRATED SWELL AND SEA DATA FILE.
      IF(fflag21)     CALL gsfile (iu06, iu21, 1, idrn, 'SWE', 'S')
!           SAVE SPECTRA FILE.
      IF(fflag25)     CALL gsfile (iu06, iu25, 1, idrn, 'SPE', 'S')
!           SAVE SWELL SPECTRA FILE.
      IF(fflag26)     CALL gsfile (iu06, iu26, 1, idrn, 'SWS', 'S')
!NEST
!           SAVE BOUNDARY VALUE FILE.
      IF( ibounc .EQ. 1 ) THEN
         iubc = 50
         DO nest = 1, nests
            IF( nest .LE. 9 ) THEN
               WRITE(cbfile, '(2x,i1)') nest
               cbfile(1:2) = 'CB'
            ELSEIF( nest .GE. 10 ) THEN
               WRITE(cbfile, '(1x,i2)') nest
               cbfile(1:2) = 'C'
            ENDIF
            IF( nest .EQ. 1 )          cbfile(1:3) = 'CBO'
            CALL gsfile (iu06, iubc, 1, idrn, cbfile, 'S')
            IF( idtpro .LT. idatee )   CALL headbc(nbounc, idelpro, fr(1), iubc, iu06, nest)
            iubc = iubc + 1
         ENDDO
      ENDIF
!NEST
!*    1.8.3 SAVE RESTART FIELDS.
!           --------------------
      WRITE(iu06,*) 'ISKA IDTOUT= ', idtout, ' IDTRES= ', idtres
      IF( irest.EQ.1 .AND. idtres.EQ.idtout .OR. idrn.EQ.'000000000000' ) THEN
! OLD CONDITION:
! IF (IREST.EQ.1.AND.IDRN.EQ.'000000000000') THEN
!*          STRESS FIELD IS DUMPED FOR RESTART.
         REWIND iu15
         WRITE(iu15) idtpro, idatewo, idawifl, idatefl, idtsou
         WRITE(iu15) u10old
         WRITE(iu15) thwold
         WRITE(iu15) usold
         WRITE(iu15) tauw
         WRITE(iu15) z0old
         WRITE(iu06,*) ' '
         WRITE(iu06,*) ' SUB. WAMODEL: WIND FIELD SAVED FOR RESTART'
         WRITE(iu06,*) ' WIND FIELD WILL BE USED UNTIL ....',' IDATEWO = ', idatewo
         WRITE(iu06,*) ' NEXT WIND FILE NAME IS ...........',' IDAWIFL = ', idawifl
         WRITE(iu06,*) ' NEXT WIND FILE WILL BE ACCESSED AT',' IDATEFL = ', idatefl
         IF( igl .EQ. 1 ) THEN
            WRITE(iu12) (((fl3(ij,k,m),ij=0,niblo),k=1,nang),m=1,nfre)
            IF( itest .GE. 2 )    WRITE(iu06,*) '   SUB. WAMODEL: START OF BLOCK OUTPUT'
         ENDIF
!*          DISPOSE RECOVERY FIELES.
         WRITE(iu06,*) 'KARIEL IDATE= ', idatee, 'COUNT = ', srcount
         rest(1:4)   = 'REST'
         rest(5:10)  = idtpro(3:8)
         rest(11:11) = '/'
#ifdef linux
! --- MAKIS FOR LINUX
         CLOSE(iu12)
         CLOSE(iu15)
#endif
         WRITE(srcmd, '(73x,i1)') srcount
         mk_dir(1:40)  = 'mkdir ../../wave/wam_in_out/swamp_nest1/'
!        mk_dir(1:35)  = 'mkdir       wam_in_out/swamp_nest1/'
         mk_dir(41:51) = rest(1:11)
         srcmd(1:42)   = 'cp ../../wave/wam_in_out/swamp_nest1/*ANAL'
!        srcmd(1:37)   = 'cp       wam_in_out/swamp_nest1/*ANAL'
         srcmd(43:77)  = ' ../../wave/wam_in_out/swamp_nest1/'
!        srcmd(38:67)  = '       wam_in_out/swamp_nest1/'
         srcmd(78:88)  = rest(1:11)
         CALL system(mk_dir)
         CALL system(srcmd)
         srcount = srcount + 1
         WRITE(iu06,*) srcmd
#ifdef linux
! --- MAKIS FOR LINUX
         OPEN(iu12, FILE = blsfile, STATUS = 'OLD', FORM = 'UNFORMATTED')
         OPEN(iu15, FILE = lawfile, STATUS = 'OLD', FORM = 'UNFORMATTED')
#endif
      ENDIF
   ENDIF
!*    BRANCHING BACK TO 1.0 FOR NEXT PROPAGATION STEP.
!   call mpi_barrier(MPI_COMM_WORLD,ierr)
ENDDO

RETURN
END SUBROUTINE wamodel
