SUBROUTINE implsch (fl3, fl, diss, ijs, ijl, ig, igl, idate)
! ----------------------------------------------------------------------
!**** *IMPLSCH* - IMPLICIT SCHEME FOR TIME INTEGRATION OF SOURCE
!****             FUNCTIONS.
!     S.D.HASSELMANN.  MPI
!     H. GUENTHER AND L. ZAMBRESKY  OPTIMIZATION PERFORMED.
!     H. GUENTHER      GKSS/ECMWF   OCTOBER 1989  NEW WIND FIELD
!                                                 INTERFACE AND
!                                                 TIME COUNTING
!     P.A.E.M. JANSSEN KNMI         AUGUST  1990  COUPLED MODEL
!     H. GUENTHER      GKSS/ECMWF   JUNE    1991  NEW SEPARATION OF
!                                                  DIAG- AND PROGNOSTIC
!                                                  PART OF SPECTRUM.
!
!*    PURPOSE.
!     --------
!       THE IMPLICIT SCHEME ENABLES THE USE OF A TIMESTEP WHICH IS
!       LARGE COMPARED WITH THE CHARACTERISTIC DYNAMIC TIME SCALE.
!       THE SCHEME IS REQUIRED FOR THE HIGH FREQUENCIES WHICH
!       RAPIDLY ADJUST TO A QUASI-EQUILIBRIUM.
!**   INTERFACE.
!     ----------
!       *CALL* *IMPLSCH (FL3, FL, IJS, IJL, IG, IGL, IDATE)*
!          *FL3*    - FREQUENCY SPECTRUM(INPUT AND OUTPUT).
!          *FL*     - DIAGONAL MATRIX OF FUNCTIONAL DERIVATIVE
!          *IJS*    - INDEX OF FIRST GRIDPOINT
!          *IJL*    - INDEX OF LAST GRIDPOINT
!          *IG*     - BLOCK NUMBER
!          *IGL*    - NUMBER OF BLOCKS
!          *IDATE* - START DATE OF SOURCE FUNCTION INTEGRATION
!     METHOD.
!     -------
!       THE SPECTRUM AT TIME (TN+1) IS COMPUTED AS
!       FN+1=FN+DELT*(SN+SN+1)/2., WHERE SN IS THE TOTAL SOURCE
!       FUNCTION AT TIME TN, SN+1=SN+(DS/DF)*DF - ONLY THE DIAGONAL
!       TERMS OF THE FUNCTIONAL MATRIX DS/DF ARE COMPUTED, THE
!       NONDIAGONAL TERMS ARE NEGLIGIBLE.
!       THE ROUTINE IS CALLED AFTER PROPAGATION FOR TIME PERIOD
!       BETWEEN TWO PROPAGATION CALLS - ARRAY FL3 CONTAINS THE
!       SPECTRUM AND FL IS USED AS AN INTERMEDIATE STORAGE FOR THE
!       DIAGONAL TERM OF THE FUNCTIONAL MATRIX.
!     EXTERNALS.
!     ---------
!       *AIRSEA*    - SURFACE LAYER STRESS AND ROUGHNESS LENGTH.
!       *CREWFN*    - CREATE A WIND FILE NAME.
!       *FEMEAN*    - COMPUTATION OF MEAN FREQUENCY AT EACH GRID POINT.
!       *INCDATE*   - UPDATE DATE TIME GROUP.
!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.
!       *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.
!       *STRESSO*   - COMPUTATION NORMALISED WAVE STRESS.
!           !!!!!!! MAKE SURE THAT SINPUT IS CALLED FIRST, STRESSO
!           !!!!!!! NEXT, AND THEN THE REST OF THE SOURCE FUNCTIONS.
!     REFERENCE.
!     ----------
!       S. HASSELMANN AND K. HASSELMANN, "A GLOBAL WAVE MODEL",
!       30/6/85 (UNPUBLISHED NOTE)
! ----------------------------------------------------------------------
use params ; use fredir ; use meanpa ; use source
use stat   ; use testo  ; use units  ; use wind
#ifndef _OPENMP
USE my_omp
#endif
IMPLICIT NONE

REAL     , INTENT(OUT):: fl3(0:niblo,nang,nfre), fl(0:niblo,nang,nfre)
REAL     , INTENT(OUT):: diss(0:niblo,nang,nfre)
INTEGER  , INTENT(IN) :: ijs, ijl, ig, igl
CHARACTER, INTENT(IN) :: idate * 12

INCLUDE'globals.h'
REAL, PARAMETER :: gzpi28 = g / 28. / zpi

LOGICAL :: newread, newfile
INTEGER :: idelwh, iareas, idnin, isbr, it1, it2, iuh, ij, j, k, l, m
REAL    :: fpmh, gtmp1, gtmp2, dfmax, delf2, fmcut, delt, umin, fac, fpm, fpm4
REAL    :: delfl(nfre)
REAL,ALLOCATABLE    :: gadiag(:), temp(:,:)
INTEGER,ALLOCATABLE :: mij(:), mfmf(:)
CHARACTER (LEN=12)  :: idatewh, idateh, idatesh
REAL    :: delt5, facb, gtemp1, gtemp2, flhab

!NIKOS These were probably remains of OLD times, i strogly believe they are not needed.
!NIKOS  They are left as comments, for historical reasons :D
!EQUIVALENCE (SL (1, 3, 1), MIJ (1) )
!EQUIVALENCE (SL (1, 5, 1), MFMF (1) )
!EQUIVALENCE (SL (1, 7, 1), GADIAG (1) )
!EQUIVALENCE (SL (1, 9, 1), TEMP (1, 1) )
#ifdef _OPENMP
INTEGER  :: omp_get_thread_num
#endif

DATA newfile / .FALSE. /

ALLOCATE(mij(niblo))    ; ALLOCATE(mfmf(niblo))
ALLOCATE(gadiag(niblo)) ; ALLOCATE(temp(niblo,nfre))

!*    1. INITIALISATION.
!        ---------------
!*    1.1 INTERNAL TIME COUNTER AND WIND FILE UNIT.
!         -----------------------------------------
idateh  = idate
idatewh = idatewo
idnin   = iu17
idatesh = idtsou
!*    1.2 WINDS ARE TAKEN FROM INTERMEDIATE STORAGE.
!         ------------------------------------------
newread = .FALSE.
IF( idateh .LT. idatewh ) THEN
!$OMP PARALLEL private(ij)
   DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
      u10new(ij) = u10old(ij,ig)
      thwnew(ij) = thwold(ij,ig)
   ENDDO
!$OMP END PARALLEL
ENDIF
! ----------------------------------------------------------------------
!*    2. COMPUTATION OF IMPLICIT INTEGRATION.
!        ------------------------------------
!         INTEGRATION IS DONE FROM IDATE UNTIL IDTPRO FOR A BLOCK
!         OF LATITUDES BETWEEN PROPAGATION CALLS.
2000 CONTINUE
!*    2.1 NEW WIND INPUT.
!         ---------------
IF( idateh .GE. idatewh ) THEN
!*    2.1.1 WINDS HAVE TO BE READ FROM FILE.
!           --------------------------------
   IF( idateh .GE. idatefl ) THEN
!*    2.1.1.1 A NEW WIND FILE HAS TO BE USED.
!             -------------------------------
      IF( ig .EQ. 1 ) THEN
!*    2.1.1.1.1 A NEW WIND FILE HAS TO BE ACCESSED.
!               -----------------------------------
         CALL crewfn (iu18, idawifl)
         newfile = .TRUE.
      ENDIF
      idnin = iu18
   ENDIF
!*    2.1.2 NEW WINDS ARE READ IN.
!           ----------------------
   READ(idnin) ida, iig, (u10new(l), l=1, niblo), (thwnew(l), l=1, niblo)
   CALL incdate (idatewh, idelwo)
   newread = .TRUE.
ENDIF
! CHANGES FOR SMALL PROP TIME STEP AND LARGE SOURCE TIME STEP
! MADE BY WEIMIN LUO, APRIL 1996
!
!* 2.1.3 WHEN IDTSOU > IDATEH, NO SOURCE TERM INTEGRATION.
!        WHEN IDTSOU <= IDATEH, SOURCE TERMS ARE INTEGRATED AND
!        UPDATE INTERNAL TIME COUNTER FOR SOURCE TERM INTEGRATION.
!
IF( idtsou .LE. idateh ) THEN
   CALL incdate (idatesh, idelt)
   IF( ig .EQ. igl ) THEN
      idtsou = idatesh
   ENDIF
   GOTO 2111
ELSE
   GOTO 3000
ENDIF
!*    2.1.3 UPDATE INTERNAL TIME COUNTER.
!           -----------------------------
2111 CONTINUE
CALL incdate (idateh, idelt)
!*    2.1.4 CALCULATE ROUGHNESS LENGTH AND FRICTION VELOCITIES.
!           ---------------------------------------------------
!$OMP PARALLEL private(m, k, ij)

!*    2.2 COMPUTE MEAN PARAMETERS.
!         ------------------------
!NIKOS CALL FEMEAN (FL3, IJS, IJL)
CALL semean (fl3, itask(1,omp_get_thread_num()+1),  &
                  itask(2,omp_get_thread_num()+1))
CALL femean (fl3, itask(1,omp_get_thread_num()+1),  &
                  itask(2,omp_get_thread_num()+1))
! ----------------------------------------------------------------------
!*    2.3 COMPUTATION OF SOURCE FUNCTIONS.
!         --------------------------------
!*    2.3.1 INITIALISE SOURCE FUNCTION AND DERIVATIVE ARRAY.
!           ------------------------------------------------
DO m = 1, nfre
   DO k = 1, nang
!NNIKOS      DO ij = 0, niblo
      DO ij = itask(1,omp_get_thread_num()+1)-1, itask(2,omp_get_thread_num()+1)
         sl(ij, k, m) = 0.
         fl(ij, k, m) = 0.
      ENDDO
   ENDDO
ENDDO
!$OMP BARRIER
!*    2.3.2 ADD SOURCE FUNCTIONS AND WAVE STRESS.
!           -------------------------------------
!NIKOS   CALL SINPUT(FL3, FL, 1, NIBLO, IG)
! iterative procedure for wave induced stress convergence
!do k = 1,5
!do k = 1,8
do k = 1,8
CALL airsea (u10new, tauw(1,ig), usnew, z0new, &
             itask(1,omp_get_thread_num()+1),  &
             itask(2,omp_get_thread_num()+1)  )
CALL sinput (fl3, fl, itask(1,omp_get_thread_num()+1),     &
                      itask(2,omp_get_thread_num()+1), ig)
CALL stresso (fl3, itask(1,omp_get_thread_num()+1),        &
                   itask(2,omp_get_thread_num()+1), ig)

end do


CALL snonlin (fl3, fl, itask(1,omp_get_thread_num()+1),    &
                       itask(2,omp_get_thread_num()+1), ig)
CALL sdissip (fl3, fl, diss, itask(1,omp_get_thread_num()+1),    &
                       itask(2,omp_get_thread_num()+1))
!SHALLOW
IF( ishallo .NE. 1 )   CALL sbottom (fl3, fl, itask(1,omp_get_thread_num()+1),    &
                                              itask(2,omp_get_thread_num()+1), ig)
! WAVES BREAKING IN SHALLOW WATERS
#ifndef waves_breaking
isbr = 0
#else
isbr = 1
IF( (ishallo .NE. 1) .AND. (isbr .EQ. 1) ) THEN
    CALL sfbrk (fl3, fl, itask(1,omp_get_thread_num()+1),    &
                         itask(2,omp_get_thread_num()+1), ig)
ENDIF
#endif
!$OMP END PARALLEL
!SHALLOW
! ----------------------------------------------------------------------
!*    2.4 COMPUTATION OF NEW SPECTRA.
!         ---------------------------
!       INCREASE OF SPECTRUM IN A TIME STEP IS LIMITED TO A FINITE
!       FRACTION OF A TYPICAL USTAR*F**(-4) SPECTRUM.
!
#ifndef newlimiter
2400 CONTINUE
fmcut = (0.1 * g**2 * 0.1e-10) / ( (fr(nfre) * zpi) **5)
delt  = idelt
fac = 3.0e-07 * g * delt * fr(nfre)
DO m = 1, nfre
   delfl(m) = fac / fr(m) **4
   umin = gzpi28 / fr(m)
   DO k = 1, nang
      DO ij = ijs, ijl
         gtmp1 = AMAX1 (1.0, 1.0 - delt * fl(ij, k, m) )
         gtmp2 = delt * sl(ij, k, m) / gtmp1
         dfmax = delfl(m) * AMAX1 (usnew(ij), umin)
         delf2 = ABS (gtmp2 + dfmax) - ABS (gtmp2 - dfmax)
         fl3(ij, k, m) = amax1 (fl3(ij, k, m) + 0.5 * delf2, fmcut)
      ENDDO
   ENDDO
ENDDO
#else
2400 CONTINUE
fmcut = (0.1 * g**2 * 0.1e-10) / ( (fr (nfre) * zpi) **5)
delt  = idelt
delt5 = 0.5 * delt
facb  = 3.0e-7 * g * delt
DO m = 1, nfre
   delfl(m) = facb * fr(m) ** ( - 4)
   DO k = 1, nang
      DO ij = ijs, ijl
         umin   = MAX (usnew(ij), gzpi28 / fmean(ij) )
         dfmax  = delfl(m) * umin * fmean(ij)
         gtemp1 = MAX ( (1. - delt5 * fl(ij, k, m) ), 1.)
         gtemp2 = delt * sl(ij, k, m) / gtemp1
         flhab  = ABS (gtemp2)
         flhab  = MIN (flhab, dfmax)
         fl3(ij, k, m) = fl3(ij, k, m) + SIGN (flhab, gtemp2)
         fl3(ij, k, m) = MAX (fl3(ij, k, m), fmcut)
      ENDDO
   ENDDO
ENDDO
#endif
! ----------------------------------------------------------------------
!*    2.5 REPLACE DIAGNOSTIC PART OF SPECTRA BY A F**(-5) TAIL.
!         -----------------------------------------------------
2500 CONTINUE
!*    2.5.1 COMPUTE MEAN PARAMETERS.
!           ------------------------
!$OMP PARALLEL private(fpmh,fpm,fpm4,ij,m,k)
CALL semean (fl3, itask(1,omp_get_thread_num()+1), &
                  itask(2,omp_get_thread_num()+1))
CALL femean (fl3, itask(1,omp_get_thread_num()+1), &
                  itask(2,omp_get_thread_num()+1))
CALL dismean (diss, itask(1,omp_get_thread_num()+1),  &
                  itask(2,omp_get_thread_num()+1))
!*    2.5.2 COMPUTE LAST FREQUENCY INDEX OF PROGNOSTIC PART OF SPECTRUM.
!*          FREQUENCIES LE MAX(4*F(PM) , 2.5*FMEAN).
!           ------------------------------------------------------------
fpmh = 2.5 / fr(1)
fpm  = 4. * gzpi28 / fr(1)
DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
!           FPM4 = FPM/(USNEW(IJ)+0.1E-90)
   fpm4     = fpm / (usnew(ij) + 0.1e-32)
   mij(ij)  = ALOG10 (fpm4) * 24.1589 + 2.
   fpm4     = fmean(ij) * fpmh
   mfmf(ij) = ALOG10 (fpm4) * 24.1589 + 1.
ENDDO
DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
   mij(ij) = MAX (mfmf(ij), mij(ij) )
   mij(ij) = MIN (mij(ij), nfre)
ENDDO
!*    2.5.3 COMPUTE TAIL ENERGY RATIOS.
!           ---------------------------
!$OMP SINGLE
DO m = 1, nfre
   delfl(m) = (1. / fr(m) ) **5.
ENDDO
!$OMP END SINGLE
DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
   gadiag(ij) = fr(mij(ij) ) **5.
ENDDO
!*    2.5.4 MERGE TAIL INTO SPECTRA.
!           ------------------------
DO m = 1, nfre
   DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
      fconst(ij,m) = 0.
      temp(ij,m)   = gadiag(ij) * delfl(m)
   ENDDO
ENDDO
DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
!NNIKOS   j = mij(ij)
!NNIKOS   DO m = 1, j
   DO m = 1, mij(ij)
      fconst(ij,m) = 1.
      temp(ij,m)   = 0.
   ENDDO
ENDDO

DO k = 1, nang
   DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
      gadiag(ij) = fl3(ij, k, mij(ij) )
   ENDDO
   DO m = 1, nfre
      DO ij = itask(1,omp_get_thread_num()+1), itask(2,omp_get_thread_num()+1)
         fl3(ij, k, m) = gadiag(ij) * temp(ij,m) + fl3(ij,k,m) * fconst(ij,m)
      ENDDO
   ENDDO
ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Consider putting this here and setting iterations=1
  CALL sinput (fl3, fl, itask(1,omp_get_thread_num()+1),     &
                        itask(2,omp_get_thread_num()+1), ig)
  CALL stresso (fl3, itask(1,omp_get_thread_num()+1),        &
                     itask(2,omp_get_thread_num()+1), ig)
  CALL airsea (u10new, tauw(1,ig), usnew, z0new, &
               itask(1,omp_get_thread_num()+1),  &
               itask(2,omp_get_thread_num()+1)  )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!$OMP END PARALLEL
! ----------------------------------------------------------------------
!*    2.6 SAVE WINDS INTO INTERMEDIATE STORAGE.
!         -------------------------------------
2600 CONTINUE
IF( newread ) THEN
   DO ij = ijs, ijl
      u10old(ij,ig) = u10new(ij)
      thwold(ij,ig) = thwnew(ij)
   ENDDO
   NEWREAD = .FALSE.
ENDIF
DO ij = ijs, ijl
   usold(ij,ig) = usnew(ij)
   z0old(ij,ig) = z0new(ij)
ENDDO
!     IF TIME LEFT BRANCH BACK TO 2.0
IF( idateh .LT. idtpro ) GOTO 2000
!     END OF TIME LOOP ALL TIME STEPS DONE
! ----------------------------------------------------------------------
!*    3. UPDATE WIND COUNTERS IF LAST BLOCK HAS BEEN DONE.
!        -------------------------------------------------
3000 IF( ig .EQ. igl ) THEN
   IF( newfile ) THEN
!*    3.1 A NEW WIND FILE HAS BEEN USED.
!         ------------------------------
!*    3.1.1 CLOSE WIND FILES SO THEY DO NOT ACCUMULATE.
!           -------------------------------------------
      CLOSE(UNIT = iu17, STATUS = 'DELETE')
!*    3.1.2 UPDATE WIND FILE TIME COUNTER AND UNITS.
!           ----------------------------------------
      iuh     = iu17
      iu17    = iu18
      iu18    = iuh
      newfile = .FALSE.
      idelwh  = MAX (idelwi, idelpro)
      CALL incdate (idawifl, idelwh)
      CALL incdate (idatefl, idelwh)
   ENDIF
!*    3.2 UPDATE WIND FIELD COUNTER.
!         --------------------------
   idatewo = idatewh
ENDIF
RETURN
END SUBROUTINE implsch
