module mod_waves_init
      !mod_waves_init.F
      !Author: Timothy Williams
      !Date:   20130612, 12:32:49 CEST
      !! subroutines for initialising waves model:
      !!  subroutine waves_setup(day1)
      !!  subroutine waves_init
#if defined (WAVES)

      use mod_xc
      implicit none

      contains

      subroutine waves_setup(day1)

      use mod_common_wavesice
      use mod_readwaves
#if defined (WAVES_ONLY)
!     use mod_forcing_waves, only: run_waves_only
      use mod_waves_only!for pivot points (ICECONS_OBS)
      use mod_hycom_nersc!for test netcdf with c,h
#endif

      implicit none
      include 'common_blocks.h'

      real, intent(in) :: day1!!1st day (1st value of dtime)
                              !!day1=1 => date=19010101)
      real  :: aday1
#if defined (WAVES_ONLY)
      !!for test output netcdf
      character(len=*), parameter ::
     &   ncobs = 'ICECONS_INPUT/test_iceobs.nc'
      real, dimension(itdm,jtdm) :: fld
      character(len=20)          :: cstr,hstr
#endif

      aday1  = abs(day1)
      if (mnproc==1) print*,'TW: day1',aday1

      if (mnproc==1) print*,'TW (mod_waves_init.F :: waves_setup):'//
     &      ' initialising waves'
      call waves_init
      if (mnproc==1) print*,'TW (mod_waves_init.F :: waves_setup):'//
     &      ' waves initialised'

      !!GET PIVOT POINTS (NEAREST-NEIGHBOURS FOR INTERPOLATION) FOR WAVE MODELS
#if defined (WAVES_WW3)
      call read_pivots_ww3
#elif defined (WAVES_WW3_ARCTIC)
      call read_pivots_ww3_arctic
#elif defined (WAVES_WAMNSEA)
      call read_pivots_wamnsea
#endif

      swh_out  = 0.
      mwp_out  = 0.
      mwd_out  = 0.

      !! read waves from external model
      !! NB even if using -DWAVES_TEST_UNIFWAVES,
      !! use mask from external model 
      call read_waves(aday1)

#if defined (WAVES_TIME_INTERP)
      mwd2  = mwd
      mwp2  = mwp
      swh2  = swh
#endif

#if defined (WAVES_ONLY)

#if defined (WAVES_ICEOBS)

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
      !!Pivot points for concentration data
#if defined (WAVES_IO_CONC2)
      !! Observed concentrations (OSISAF),
      !! Get pivot points (for nearest neighbour interpolation)
      call read_pivots_osisaf
      cstr  = 'Conc (OSISAF)'
#elif defined (WAVES_IO_CONC3)
      !! Observed concentrations (AMSR2 3.125km),
      !! Get pivot points (for nearest neighbour interpolation)
      call read_pivots_amsr2_3125
      cstr  = 'Conc (AMSR2 3.125km)'
#else
      !! Observed concentrations (AMSR-E Bremen),
      !! Get pivot points (for nearest neighbour interpolation)
      call read_pivots_amsre6250
      cstr  = 'Conc (AMSR-E)'
#endif
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   
      !!Pivot points for thickness data, if necessary;
#if defined (WAVES_IO_THICK1)
      !!  observed thicknesses (SMOS)
      !! Get pivot points (for nearest neighbour interpolation)
      if (mnproc==1) print*,'TW: getting SMOS pivots'
      call read_pivots_smos
      hstr  = 'Thickness (SMOS)'
#else
      !!  thickness estimated (not read in)
      hstr  = 'Thickness (Estimate)'
#endif
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   

#else/* ! WAVES_ICEOBS */
      !!(else: Get c,h from daily average files)
      !!(no pivot pts needed since on same grid)
      cstr  = 'Conc (DA)'
      hstr  = 'Thickness (DA)'
#endif/* ! WAVES_ICEOBS */

#if defined (WAVES_ICEOBS)
      if (.false.) then
         !!check data is read correctly;
         call read_ice_thick_conc(aday1)
         if ((itest.ge.0.0).and.(jtest.ge.0.0)) then
            print*,'TW: i/j',itest+i0,jtest+j0
#if defined (WAVES_IO_THICK1)
            print*,'TW: i/j-piv (SMOS)',
     &         ipiv_thick(itest,jtest),jpiv_thick(itest,jtest)
#endif
            print*,'TW: lon/lat',plon(itest,jtest),plat(itest,jtest)
            print*,'TW: c/h',ficem(itest,jtest),hicem(itest,jtest)
         end if
         !!
         call xcaget(fld,plon,0)
         call ncdraft(trim(ncobs),fld,'lon','clobber')
         call xcaget(fld,plat,0)
         call ncdraft(ncobs,fld,'lat','')
         call xcaget(fld,ficem,0)
         call ncdraft(ncobs,fld,trim(cstr),'')
         call xcaget(fld,hicem,0)
         call ncdraft(ncobs,fld,trim(hstr),'')
         if (mnproc==1) print*,
     &      'TW (mod_waves_init.F): dumping observed conc/thickness 
     &   and exiting'
         call xcstop('test observed conc/thickness done')
      end if
#endif/* WAVES_ICEOBS*/

!!!!      if (mnproc==1) print*,'tw (mod_waves_init.f): running waves-only'
!!!!      call run_waves_only
!!!!      if (mnproc==1) print*,
!!!!     &   'tw (mod_restart.f): waves-only run finished'
!!!!      call xcstop('tw (mod_restart.f): waves-only run finished')
!!!!      stop
#endif/* WAVES_ONLY */

      end subroutine waves_setup

      subroutine waves_init
      !! may use subroutine 'rotate2' to get uwave,vwave here
      !! - would then need mod_xc.F and mod_hycom_nersc.F
      !use mod_xc
      !use mod_hycom_nersc
      use mod_xc
      use mod_common_wavesice
      use mod_wim_prams
      use mod_forcing_waves,only:
     &   waves_set_icemask,waves_fsd_reset

      implicit none

      include 'common_blocks.h'

      real, dimension (itdm,jtdm) :: dep_all,dx_all
      real, parameter             :: dx_large = 1.0e6!! 1000km - this number just needs
                                                     !! to be larger than all the grid cell sizes
                                                     !! in all model setups
                                                     !! (eg TP4a0.12,BS*,FR* etc);
      real   :: stress_c

      !!wave directions:
      real           :: th1,th2,dtheta
      integer        :: wth,w
      real           :: tmin,tmax,ommin,ommax,om,dw

      integer,parameter  :: ndir  = n_wavdir
      integer,parameter  :: nw    = n_wavefreq

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! set up directions first

      !!set wave directions
      th1      = avgmwd-dth_tot/2.0
      th2      = avgmwd+dth_tot/2.0
      dtheta   = (th2-th1)/ndir

      !!set integral weights wrt directions
      if (ndir.eq.1) then
         wt_theta = 1.
      else
         wt_theta = dtheta
      end if

      if(mnproc==1)print*,' '
      if(mnproc==1)print*,'--------------------------------------------'
      if(mnproc==1)print*,'wave directions'
      do wth = ndir,1,-1
         wavdir(wth) = th1+(wth-0.5)*dtheta
         if(mnproc==1)print*,wavdir(wth)
      end do
      if(mnproc==1)print*,'--------------------------------------------'
      if(mnproc==1)print*,' '

      if (ACoption.eq.0) then
         tmin  = 6.0
         tmax  = 15.9
      else
         tmin  = 1/0.4
         tmax  = 1/0.042
      end if
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! set up frequencies

      !!set wave frequencies
      ommin = 2*pi/tmax
      ommax = 2*pi/tmin
      dw    = (ommax-ommin)/(nw-1.)

      if(mnproc==1)print*,' '
      if(mnproc==1)print*,'--------------------------------------------'
      if(mnproc==1)print*,'wave periods'
      do w = 1,nw
         om             = ommin+(w-1)*dw
         wave_period(w) = 2.*pi/om
         if(mnproc==1)print*,wave_period(w)
      end do
      if(mnproc==1)print*,'--------------------------------------------'
      if(mnproc==1)print*,' '

      !! weights for integral wrt omega (Simpson's rule)
      wt_omega     = 2.0
      wt_omega(1)  = 1.0
      wt_omega(nw) = 1.0
      do w=2,nw-1,2
         wt_omega(w)  = 4.0
      end do
      wt_omega = (dw/3.0)*wt_omega
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



      stress_c = 0.0

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!Initialise ice mask & dfloe
      call waves_set_icemask

      if (FSD_INIT) then
         !if FSD not initialised in read_restart_ice_ab
         call waves_fsd_reset
      end if
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! define LAND_MASK
      !! & find min grid size (can't use minval as scpx is 0/1/? on land)
      LAND_MASK   = 1.0
      where (depths(1-nbdy:ii+nbdy,1-nbdy:jj+nbdy).gt.0.0)
     &         LAND_MASK   = 0.0

      call xcaget(dx_all,scpx,0)
      call xcaget(dep_all,depths,0)
      where (dep_all.eq.0.0) dx_all = dx_large
      scpx_min = minval(dx_all)
      if(mnproc==1)print*,'TW: dx_large, dxmin',dx_large,scpx_min
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      DO_WAVES_INIT  = .true.
#if !defined (WAVES_NOSAVE)
      sdf_out  = 0.0
      swh_out  = 0.0
      mwp_out  = 0.0
      mwd_out  = 0.0
#endif

      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !! set some parameters for the WIM
      sigc  = sigc0*exp(-5.88*sqrt(brine_vol))  !Timco & O'Brien (1994)

      !!set Young's modulus
      if (YOUNG_OPT.eq.1) then
         young = young0*(1-3.51*brine_vol)-1.0e9   !Williams et al (2013a)
      else
         young = 2.0e9
      end if

      !!set breaking strain
      if (BRK_OPT.eq.1) then
         epsc  = sigc/young
      else
         stress_c = 2.6*sigc!!=E/(1-nu^2)*strain_c = thin plate (plane stress: \sigma_33=0)
         epsc     = (1-poisson**2)*stress_c/young
      end if

!     !!parameters for Robinson-Palmer dispersion relation
!     prams_RP(1) = young
!     prams_RP(2) = poisson
!     prams_RP(3) = rhowtr
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      end subroutine waves_init
#endif/* WAVES */

      end module mod_waves_init

!     -----------------------------------------------------------------------------------
!     REVISION HISTORY
!     -----------------------------------------------------------------------------------
!
!     20151108 (rev 400):
!        Set wave frequencies here now