SUBROUTINE wavemdl
! ----------------------------------------------------------------------
!**** *WAVEMDL* - SUPERVISES EXECUTION OF THE WAVE MODEL.
!     LIANA ZAMBRESKY    GKSS/ECMWF    OCTOBER 1988
!     MODIFIED BY H. GUNTHER   ECMWF   MARCH 1990
!     P. KATSAFADOS       HUA/HCMR     JULY 2014 SENDS SWH, CHARNOCK etc TO THE ATMOS MODEL
!*    PURPOSE.
!     --------
!       THIS SUB SUPERVISES THE EXECUTION OF THE
!       MAIN MODULES FOR WAM MODEL INITIALIZATION,
!       WIND FIELD PREPROCESSING,
!       WAM MODEL EXECUTION.
!**   INTERFACE.
!     ----------
!       SEE MAIN MODULES SUB INITMDL, PREWIND AND WAMODEL.
!     METHOD.
!     -------
!       THE FIRST TIME WAVEMDL IS CALLED, THE WAM MODEL IS INITIALIZED.
!       THIS INITIALIZATION INCLUDES GETTING THE INITIAL SEA STATE
!       FILES, FILLING COMMON BLOCKS, DEFINING THE GRID AND SETTING
!       GENERAL PARAMETERS.
!       IN THE FIRST AND ALL SUBSEQUENT CALLS TO WAVEMDL PREWIND
!       REFORMATS THE WINDS INTO THE WAM MODEL BLOCKED STRUCTURE AND
!       THE WAM MODEL IS EXECUTED.
!       EACH CALL TO WAMODEL INTEGRATES THE WAVE SPECTRA FORWARD
!       IN TIME BY ONE INPUT WIND TIME STEP OR PROPAGATION TIME
!       STEP, WHAT EVER IS GREATER.
!     EXTERNALS.
!     ----------
!       *INITMDL*   - INITIALIZES THE WAM MODEL.
!                     GETS RECOVERY FILES.
!                     SETS COMMON BLOCKS NECESSARY TO DEFINE
!                     THE GRID AND BLOCKING STRUCTURE.
!                     DEFINES GENERAL PARAMETERS.
!       *PREWIND*   - REFORMATS WINDS ON THE GAUSSIAN GRID
!                     INTO THE WAM MODEL BLOCKED STRUCTURE.
!       *WAMODEL*   - INTEGRATES THE WAVE SPECTRA FORWARD IN TIME BY
!                     ONE WIND INPUT TIME STEP OR ONE PROPAGATION
!                     TIME STEP, WHATEVER IS GREATER.
!     REFERENCE.
!     ----------
!       NONE.
! ----------------------------------------------------------------------
USE params ; USE wind1 ; USE testo ; USE ubuf  ; use wind ; use map
USE stat ; USE meanpa
#ifdef assimilation
USE stat
#endif
use couple_mpi
use shallow ; use spe ; use coupl
IMPLICIT NONE


INCLUDE 'globals.h'
REAL    :: betag(ngx,ngy),z0g(ngx,ngy),swhg(ngx,ngy)
REAL    :: cpp(niblo), om, age, ageg(ngx,ngy), wdir(ngx,ngy), wmask(ngx,ngy)
REAL    :: ustar(ngx,ngy)
CHARACTER (LEN=12),SAVE  :: IDTWIR
CHARACTER (LEN=14) :: OUTBETANAME 
#ifdef assimilation
CHARACTER (LEN=12) :: assim_dates(100), assim_date  
#endif
LOGICAL :: frstime  , frsbeta , liu_theory
!     NADV  NUMBER OF PROPAGATION TIME STEPS IN ONE CALL TO WAMODEL.
INTEGER :: nadv, nsubadv
INTEGER :: k, ij, nn, ijs, ijl, ijss, ijsst, ijllt, ijll, nass, jtask
INTEGER :: ix, iy
INTEGER :: ii, jj
REAL    :: aminu, amaxu, aminw, amaxw, aminus, amaxus, amind, amaxd
! ----------------------------------------------------------------------
DATA frstime, nadv, nsubadv / .TRUE., 0, 1 /
!*    1. THE FIRST CALL TO WAVEMDL PERFORMS INITIALIZATION.
!        --------------------------------------------------

!     CALL MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)

   frsbeta = .false.

! ----- Coupling theory according to Liu et al. 2011 paper. Activated
!       when liu_theory set to true
   liu_theory = .false.


IF(frstime) THEN  
   frsbeta = .true.
   CALL initmdl(nadv, nsubadv)
! ---- FL1 comes in from model initialization
   CALL wavenumber(fl1, cpp)
   IDTWIR = IDATEA
   betag=9999.
   swhg=0.
   wmask = 0.
   ustar = 0.
   DO ij = 1, niblo
   IX = IXLG (IJ, 1)
! write from north to south
!  IY = NGY - KXLT (IJ, 1) + 1
! write from  south to north
   IY =  KXLT (IJ, 1) + 1
   betag(ix,iy)    = G*Z0OLD(IJ,1)/ (usold(ij,1) **2 + 0.0001)
   z0g(ix,iy)    = Z0OLD(IJ,1)
   swhg(ix,iy)    = 4.* SQRT(MAX(EMEAN(IJ),0.))
   wdir(ix,iy)    = DEG*THQ(IJ)
   wmask(ix,iy)    = 1.
   ageg(ix,iy) = CPP(IJ)/(USOLD(IJ,1)+0.00001)
   ustar(ix,iy) = usold(ij,1)

   if (liu_theory) then
   om = min(1., 0.64/(0.41*(USOLD(IJ,1)+0.00001)))
   age = CPP(IJ)/(USOLD(IJ,1)+0.00001)
   print*,'om, age, usold,xkappa init: ',om,age,USOLD(IJ,1),xkappa
   if (age.lt. 35.) then
    betag(ix,iy) = (0.085*age**1.5)**(1.-1./om)* \
                   (0.03*age*exp(-0.14*age))**(1./om)
   else
    betag(ix,iy) = 17.60**(1.-1./om)*0.008**(1./om)
   end if
   end if

   END DO
! write out beta (based on z0old , usold) for the first time
!    if (myrank .EQ. 1) then
     !print *, 'WAM WAVEMDL MYPE,NPES ', mype2,npes2
      aminu=100000.
      amaxu=-100000.
      aminw=100000.
      amaxw=-100000.
      amind=100000.
      amaxd=-100000.
      do ii=1,ngx
      do jj=1,ngy
         if (betag(ii,jj).ne.9999.) then
            if (betag(ii,jj).lt.aminu) aminu=betag(ii,jj)
            if (betag(ii,jj).gt.amaxu) amaxu=betag(ii,jj)
         endif
         if (swhg(ii,jj).gt.1.E-5 .and. swhg(ii,jj).lt.1.E+1) then
            if (swhg(ii,jj).lt.aminw) aminw=swhg(ii,jj)
            if (swhg(ii,jj).gt.amaxw) amaxw=swhg(ii,jj)
         endif
!         if (wdir(ii,jj).gt.1.E-5) then
            if (wdir(ii,jj).lt.amind) amind=wdir(ii,jj)
            if (wdir(ii,jj).gt.amaxd) amaxd=wdir(ii,jj)
!         endif
      enddo
      enddo
      write(6,*)'WAM SENDS INITIAL betag min-max ',aminu,amaxu
      write(6,*)'WAM SENDS INITIAL swhg min-max ',aminw,amaxw
      write(6,*)'WAM SENDS INITIAL wdir min-max ',amind,amaxd
!
     call mpi_send(idtpro,12,MPI_CHARACTER,atmos_proc_c,tagi, \
                   MPI_COMM_WORLD,status,ierr)
     call mpi_send(betag,ngx*ngy,MPI_REAL,atmos_proc_c,tagw, \
                   MPI_COMM_WORLD,status,ierr)
! PETROS Send  z0 and Significant wave height
    call mpi_send(z0g,ngx*ngy,MPI_REAL,atmos_proc_c,tagz, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(swhg,ngx*ngy,MPI_REAL,atmos_proc_c,tagh, \
                  MPI_COMM_WORLD,status,ierr)
! New variables to be sent (wdir, mask, wave age)
    call mpi_send(wdir,ngx*ngy,MPI_REAL,atmos_proc_c,tagd, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(wmask,ngx*ngy,MPI_REAL,atmos_proc_c,tagm, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(ageg,ngx*ngy,MPI_REAL,atmos_proc_c,taga, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(ustar,ngx*ngy,MPI_REAL,atmos_proc_c,tags, \
                  MPI_COMM_WORLD,status,ierr)
!     CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
     print *,'WAM WAVEMDL SENDS INITDATE,betag,z0g,swhg(ngx/2,ngy/2)', \
            IDTPRO,betag(ngx/2,ngy/2),z0g(ngx/2,ngy/2),swhg(ngx/2,ngy/2)
!    endif


   frstime = .FALSE.  
   IF( itest .GE. 1 )  WRITE(IU06,*) ' SUB. WAVEMDL: INITMDL DONE'  
#ifdef assimilation
   nass = 0
   OPEN(66, FILE='ASSIMILATION_DATES', STATUS='OLD', ERR=7000)
   READ(66,*) nass
   DO k = 1, nass
      READ(66, '(A12)') assim_dates(k)
   ENDDO
   CLOSE(66)
   7000 CONTINUE
#endif
ENDIF

! SEND TO ATMOSPHERIC MODEL BETA DATA                       
if (idtpro .eq. idtwir) then
  if (.not. frsbeta) then
   CALL wavenumber(fl3, cpp)
   betag=9999.
   swhg=0.
   wmask=0.
   ustar=0.
   aminus=100000.
   amaxus=-100000.
   DO ij = 1, niblo
   IX = IXLG (IJ, 1)
! write from north to south
!  IY = NGY - KXLT (IJ, 1) + 1
! write from  south to north
   IY =  KXLT (IJ, 1) + 1
!  betag(ix,iy)    = G*Z0NEW(IJ)/ (usnew(ij) **2 + 0.0001)
   betag(ix,iy)    = G*Z0NEW(IJ)/ max(usnew(ij) **2 , 1.e-3)
   z0g(ix,iy)    = Z0NEW(IJ)
!   if (emean(ij).lt.0.) emean(ij)=0.
   swhg(ix,iy)    = 4.* SQRT(MAX(EMEAN(IJ),0.))
   wdir(ix,iy)    = DEG*THQ(IJ)
   wmask(ix,iy)    = 1.
   ageg(ix,iy) = CPP(IJ)/(USNEW(IJ)+0.00001)
   ustar(ix,iy) = usnew(ij)
   if (USNEW(ij).lt.aminus) aminus=USNEW(ij)
   if (USNEW(ij).gt.amaxus) amaxus=USNEW(ij)

   if (liu_theory) then
   om = min(1., 0.64/(0.41*(USNEW(IJ)+0.00001)))
   age = CPP(IJ)/(USNEW(IJ)+0.00001)
   if (age.lt. 35.) then
    betag(ix,iy) = (0.085*age**1.5)**(1.-1./om)* \
                   (0.03*age*exp(-0.14*age))**(1./om)
   else
    betag(ix,iy) = 17.60**(1.-1./om)*0.008**(1./om)
   end if
   end if

   END DO
   print *, 'WAM min-max USNEW: ', aminus,amaxus

!    if (myrank .EQ. 1) then
     print *, 'WAM WAVEMDL SENDS DATE: ', IDTPRO
     aminu=100000.
     amaxu=-100000.
     aminw=100000.
     amaxw=-100000.
     amind=100000.
     amaxd=-100000.
     do ii=1,ngx
     do jj=1,ngy
        if (betag(ii,jj).ne.9999.) then
           if (betag(ii,jj).lt.aminu) aminu=betag(ii,jj)
           if (betag(ii,jj).gt.amaxu) amaxu=betag(ii,jj)
        endif
        if (swhg(ii,jj).gt.1.E-5 .and. swhg(ii,jj).lt.1.E+1) then
           if (swhg(ii,jj).lt.aminw) aminw=swhg(ii,jj)
           if (swhg(ii,jj).gt.amaxw) amaxw=swhg(ii,jj)
        endif
!         if (wdir(ii,jj).gt.1.E-5) then
            if (wdir(ii,jj).lt.amind) amind=wdir(ii,jj)
            if (wdir(ii,jj).gt.amaxd) amaxd=wdir(ii,jj)
!         endif
!        print*,ii,jj,betag(ii,jj),swhg(ii,jj),wdir(ii,jj),wmask(ii,jj)
!        print*,ii,jj,swhg(ii,jj),wmask(ii,jj),ageg(ii,jj)
     enddo
     enddo
     write(6,*)'WAM SENDS betag min-max ',aminu,amaxu
     write(6,*)'WAM SENDS swhg min-max ',aminw,amaxw
     write(6,*)'WAM SENDS wdir min-max ',amind,amaxd
     call mpi_send(idtpro,12,MPI_CHARACTER,atmos_proc_c,tagi, \
                   MPI_COMM_WORLD,status,ierr)
     call mpi_send(betag,ngx*ngy,MPI_REAL,atmos_proc_c,tagw, \
                   MPI_COMM_WORLD,status,ierr)
! PETROS Send  z0 and Significant wave height
    call mpi_send(z0g,ngx*ngy,MPI_REAL,atmos_proc_c,tagz, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(swhg,ngx*ngy,MPI_REAL,atmos_proc_c,tagh, \
                  MPI_COMM_WORLD,status,ierr)
! New variables to be sent (wdir, mask, wave age)
    call mpi_send(wdir,ngx*ngy,MPI_REAL,atmos_proc_c,tagd, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(wmask,ngx*ngy,MPI_REAL,atmos_proc_c,tagm, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(ageg,ngx*ngy,MPI_REAL,atmos_proc_c,taga, \
                  MPI_COMM_WORLD,status,ierr)
    call mpi_send(ustar,ngx*ngy,MPI_REAL,atmos_proc_c,tags, \
                  MPI_COMM_WORLD,status,ierr)

!     CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
     print *, 'WAM WAVEMDL SENDS DATE,betag,z0g,swhg(ngx/2,ngy/2): ', \
            IDTPRO,betag(ngx/2,ngy/2),z0g(ngx/2,ngy/2),swhg(ngx/2,ngy/2)
!    endif
  end if
    call incdate(idtwir,idelwi)
endif
! ----------------------------------------------------------------------
!*    2. INTEGRATE THE WAVE SPECTRA FORWARD IN TIME.
!        -------------------------------------------
!*    2.1 REFORMAT WINDS FROM GAUSSIAN TO BLOCKED.
!         ----------------------------------------
CALL prewind
IF( itest .GE. 1 )    WRITE(iu06,*) ' SUB. WAVEMDL: PREWIND DONE'  
!*    2.2 INTEGRATE THE WAVE SPECTRA FORWARD IN TIME.
!         -------------------------------------------
CALL wamodel(nadv,nsubadv)
IF( itest .GE. 1 )    WRITE(iu06,*) ' SUB. WAVEMDL: WAMODEL DONE'


RETURN
END SUBROUTINE wavemdl
