!==============================================================================

MODULE W3GDATMD  ! reduced version  of W3GDATMD.FTN

PUBLIC

   INTEGER,DIMENSION(:,:),ALLOCATABLE            :: TRIGP
   DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE   :: XYB
   INTEGER                                       :: NTRI, NX, NK
   REAL, ALLOCATABLE                             :: FREQ(:), FREQ1(:), FREQ2(:), DF(:)
      INTEGER, PARAMETER      :: RLGTYPE = 1
      INTEGER, PARAMETER      :: CLGTYPE = 2
      INTEGER, PARAMETER      :: UNGTYPE = 3   

END MODULE W3GDATMD

!==============================================================================

!=================================================
!
	PROGRAM ATMOSPHERIC_LP12_RESPONSE
!
!               +----------------------------------------+
!               | Equivalent of matlab code in fortran   |
!               +----------------------------------------+
!==================================================

! INPUTS PARAMETERS (in file .inp)
! * name of stations (or number of stations)
! * type of propa (for now just 1)
! * Name of paths : - PathNC for wave files
!		    - SavePath
!		    - WindPath
! 		    - DistazPath (for summation of cells)
! 		    - PowerPath (and name)
! * STD Name of files
! * AZI boundaries

!
  USE C_CALDAT
  USE W3SERVMD
  USE W3GDATMD
  USE W3TRIAMD
  USE W3GSRUMD
  USE NETCDF
  USE ieee_arithmetic
!
  IMPLICIT NONE
!
  INTEGER                          :: NDSE, NDSO, NDSI
!
  CHARACTER*1                      :: COMSTR
  CHARACTER*4                      :: CSTRG
  REAL, PARAMETER                  :: g    =9.81     ! gravity          in m/s^2
  REAL, PARAMETER                  :: rhow =1026     ! density of water in kg/m^3
  REAL, PARAMETER                  :: rhos =2600     ! density of rock  in kg/m^3
  REAL, PARAMETER                  :: betas=2800     ! shear wave velocity of rock in m/s
  REAL, PARAMETER                  :: GRAV = 9.81
  REAL, PARAMETER                  :: pi=3.1415926536, R_E=4E7/(2*pi) ! Earth radius in m
  DOUBLE PRECISION, PARAMETER      :: lg10=alog(10.)
! General variables (input reading)
!.........................................................
  CHARACTER*128			   :: varp2l, REFTAG, WINDTAG
  CHARACTER*256                    :: filein, pathp2l, filepower, filepowerbis, filewind, filestagen
  CHARACTER(LEN=10), ALLOCATABLE   :: STATIONS(:)
  INTEGER                          :: date1(6), date2(6), dimid(4), dimln(4), NTH, TH1, id, mm, iyyy
  INTEGER(KIND=4)                  :: JDAY1, JDAY2, JDAYREF
  INTEGER                          :: IRET, NCID, VARIDTMP, NCIW
  INTEGER                          :: I
  INTEGER                          :: NSTATIONS, NLONMOD, NLATMOD, NFMOD, NT, NWINDAZI, NTIMEWIND
  LOGICAL                          :: FLAG_POWBIS
  DOUBLE PRECISION                 :: DT, DTMOD ! Time step
  REAL, ALLOCATABLE                :: LONO(:), LATO(:) ! LON and LAT of stations
  REAL, ALLOCATABLE                :: MODLAT(:), MODLON(:), AZIWIND(:), MODFREQ(:)
  REAL, ALLOCATABLE                :: VRATIO(:,:)
  REAL, POINTER                    :: ALON(:,:),ALAT(:,:)
  REAL                             :: latmin, latmax, lonmin, lonmax
  DOUBLE PRECISION, ALLOCATABLE    :: PMAP(:,:,:), PMAPBIS(:,:,:)
  DOUBLE PRECISION, ALLOCATABLE    :: TIMES(:), TIMEWIND(:)
  DOUBLE PRECISION                 :: WTMIN, WTMAX, DTW, MAXDISTANCE
! Variables for station file handling
!...............................................................
  CHARACTER*256                    :: filesta
  INTEGER                          :: NCIS
  INTEGER                          :: NSTAAZI, NILON, NILAT, NWEIGHT, NDIS, NSTARTAZI, NENDAZI
  INTEGER, ALLOCATABLE             :: ILON(:),ILAT(:), STARTAZI(:), ENDAZI(:), ILON0(:), ILAT0(:)
  DOUBLE PRECISION                 :: DTH
  REAL, ALLOCATABLE                :: AZISTA(:)
  DOUBLE PRECISION, ALLOCATABLE    :: WEIGHT(:)
! Variables for source file handling
!...............................................................
  CHARACTER*256                    :: filenc, filencold
  CHARACTER*128                    :: UNITP2L, UNITP2LTMP
  INTEGER                          :: NCIP, DIMT, DIMF, VARF, VARTIME, VARIDP2L
  INTEGER                          :: NTMOD, NLONNC,NLATNC, LENU, NTRMAX,NTI, NTR, STEPMOD, L1, L2, IBEGINMOD, pp, iin
  DOUBLE PRECISION, allocatable    :: indmod(:) ! before it was integer
  INTEGER                          :: ITR, ITIME, IMIN, IMAX,indfreqdeb
  REAL                             :: SCALEP2L, FILLP2L, SCALEP2LTMP, FILLP2LTMP
  REAL, ALLOCATABLE                :: MODFREQB(:),FREQB(:)!FREQ(:), FREQ1(:), FREQ2(:), DF(:)
  DOUBLE PRECISION                 :: NOW, NOWNEW
  DOUBLE PRECISION, ALLOCATABLE    :: TIMESMOD(:),DISTANCE(:), TABP2F(:)
  DOUBLE PRECISION, ALLOCATABLE    :: P2FT(:,:,:,:)
! Variables in the do-loop
!...............................................................
  INTEGER                          :: IT, K, AZ, STARTA, ENDSTA, LL
  INTEGER                          :: PTDIR, IO, IA, MAXIAMP,AZS2
  REAL                             :: AZI, AZE, AZS, VRAZE
  REAL, ALLOCATABLE                :: VRATIOAZ(:)
  LOGICAL, ALLOCATABLE             :: SOURCEOK(:)!,MASK_DIST(:)
  DOUBLE PRECISION, ALLOCATABLE    :: WEIGHT2(:), DIS(:),DIS_MASKED(:), SOURCE(:),SOURCEBIS(:)
  DOUBLE PRECISION, ALLOCATABLE    :: ATTTMP1(:), AMP1(:), ATTTMP2(:), AMP2(:)
  INTEGER                          :: NKK,indf1,indf2
  REAL                             :: limf1,limf2
  real, allocatable                :: dfind(:)
  DOUBLE PRECISION, ALLOCATABLE    :: F_spec1(:,:,:), F_spec2(:,:,:), F_specBIS1(:,:,:), F_specBIS2(:,:,:)
  Double precision, allocatable    :: Amp_t_az_Delf1(:,:,:),Amp_t_az_Delf2(:,:,:),Amp_t_az_DelfBIS1(:,:,:),Amp_t_az_DelfBIS2(:,:,:)
  Double precision, allocatable    :: f_t_az_mamp1(:,:),f_t_az_mamp2(:,:),f_t_az_mampBIS1(:,:),f_t_az_mampBIS2(:,:)
  Double precision, allocatable    :: Az_t_f_mamp1(:,:),Az_t_f_mamp2(:,:),Az_t_f_mampBIS1(:,:),Az_t_f_mampBIS2(:,:)
  Double precision, allocatable    :: Amp_t_az_sum1(:,:),Amp_t_az_sum2(:,:),Amp_t_az_sumBIS1(:,:),Amp_t_az_sumBIS2(:,:)
  Double precision, allocatable    :: Amp_t_az_mf1(:,:),Amp_t_az_mf2(:,:),Amp_t_az_mfBIS1(:,:),Amp_t_az_mfBIS2(:,:)
  Double precision, allocatable    :: Amp_t_f_sum1(:,:),Amp_t_f_sum2(:,:),Amp_t_f_sumBIS1(:,:),Amp_t_f_sumBIS2(:,:)
  Double precision, allocatable    :: Amp_t_f_maz1(:,:),Amp_t_f_maz2(:,:),Amp_t_f_mazBIS1(:,:),Amp_t_f_mazBIS2(:,:)
!  DOUBLE PRECISION, ALLOCATABLE    :: F_spec_tot1(:,:), F_spec_totBIS1(:,:)
!  DOUBLE PRECISION, ALLOCATABLE    :: F_spec_tot2(:,:), F_spec_totBIS2(:,:)
!  DOUBLE PRECISION, ALLOCATABLE    :: DISTM1(:,:,:), AMPM1(:,:,:), F_spec1(:,:,:)
!  DOUBLE PRECISION, ALLOCATABLE    :: DISTM2(:,:,:), AMPM2(:,:,:), F_spec2(:,:,:)
!  DOUBLE PRECISION, ALLOCATABLE    :: DISTMBIS1(:,:,:), AMPMBIS1(:,:,:), F_specBIS1(:,:,:)
!  DOUBLE PRECISION, ALLOCATABLE    :: DISTMBIS2(:,:,:), AMPMBIS2(:,:,:), F_specBIS2(:,:,:)
!
! Saving file
!...............................................................
  CHARACTER*256                    :: filencout1, filencout2,pathNCOUT
  INTEGER                          :: ncio1, ncio2,varido1(28), varido2(28),varid(28),dimlno(6),dimido(6),dimln1(6), dimid1(6)
interface 
SUBROUTINE attenuation(DIS,f,Vratio,ATTTMP,LL)
   INTEGER                         :: LL, f1, F2
   REAL, DIMENSION(LL)             :: Vratio, beta_val
   DOUBLE PRECISION, DIMENSION(LL) :: DIS
   REAL                            :: f_table(6), Veff_table(12), alpha_table(6),sigma_table(6), beta_val0(12)
   REAL                            :: beta_table(12,6)
   REAL, DIMENSION(LL,12)          :: Vratio2,Veff_table2
   INTEGER, DIMENSION(LL)          :: V1,V2
   DOUBLE PRECISION, ALLOCATABLE   :: ATTTMP(:)
   REAL                            :: f, delta, MinVeff, MaxVeff, alpha_val,sigma_val
end subroutine attenuation

SUBROUTINE attenuationAZE (DIS,f,Vraze,ATTTMP,LL)
   INTEGER                         :: LL, f1, F2, V1,V2
   DOUBLE PRECISION, DIMENSION(LL) :: DIS
   REAL                            :: f_table(6), Veff_table(12), alpha_table(6),sigma_table(6), beta_val0(12),Vraze
   REAL                            :: beta_table(12,6)
   DOUBLE PRECISION, ALLOCATABLE   :: ATTTMP(:)
   REAL                            :: f, delta, MinVeff, MaxVeff, alpha_val, beta_val, sigma_val
end subroutine attenuationAZE

SUBROUTINE attenuationNoWind(DIS,f,ATTTMP,LL)
   INTEGER                         :: LL, f1, F2
   DOUBLE PRECISION, DIMENSION(LL) :: DIS
   REAL                            :: f_table(6), alpha_table(6),sigma_table(6)
   REAL                            :: beta_table(6)
   DOUBLE PRECISION, ALLOCATABLE   :: ATTTMP(:)
   REAL                            :: f, delta, alpha_val,sigma_val, beta_val
end subroutine attenuationNoWind

SUBROUTINE DISTAZ (la1,lo1,la2,lo2,az2from1,az1from2)
! Given 2 points : 1(la1,lo1) and 2(la2,lo2) we compute two angles: az2from1 : azimuth of point 2 seen from point 1, az1from2: azimuth of point 1 seen from point2
implicit none
REAL, INTENT(IN)          :: la1(1),lo1(1),la2,lo2
REAL, INTENT(OUT)         :: az2from1, az1from2
REAL                      :: sa,ssa,sca,so,ea,eo,sea,eca,ecs,ses,rad
end subroutine distaz
end interface
!=============================================================================== 
! 1.a IO set-up
!=============================================================================== 
  NDSI=10
  NDSE=6
  NDSO=6
  JDAYREF=julday(1,1,1990)
!===============================================================================  
! 1.b Read input parameters from files
!===============================================================================  
  WRITE(6,*) 'IARG:',iargc()
! -------- Read name of input file : given in command line or default 
  IF (iargc().GE.1) THEN
    call getarg(1,filein)
  ELSE 
    filein='atmospheric_response_example.inp'
  ENDIF
!
  WRITE (NDSO,*) 'Uses input file:',TRIM(filein)
  open(NDSI,file=filein) 
  READ (NDSI,'(A)') COMSTR
  IF (COMSTR.EQ.' ') COMSTR = '$'
! ---------------------------------------
! 1.b.1.GLOBAL INPUTS
! ---------------------------------------
! ------- Start date------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  read(NDSI,'(I4I2I2,X,3I2)') date1
  WRITE(6,*) 'date 1',date1
  JDAY1 = julday(date1(3),date1(2),date1(1)) - JDAYREF
! ------- End date------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  read(NDSI,'(I4I2I2,X,3I2)') date2
  WRITE(6,*) 'date 1', date2
  JDAY2 = julday(date2(3),date2(2),date2(1)) -JDAYREF
! ------- Time Step in days------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,*) DT
  WRITE (NDSO,*) 'Time step:',DT
! ------- Path name for P2L file------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') pathp2l
  WRITE (NDSO,*) 'Will use WW3 files starting with :',trim(pathp2l)
! ------- Reflection Tag (if needed)------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') REFTAG
! ------- Wind Tag (How to compute attenuation) ! TO BE CHANGED
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') WINDTAG
! ------- Name of NetCDF power file------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') filepower
  WRITE (NDSO,*) 'Power coefficient taken from file:',trim(filepower)
! ------- Name of NetCDF power file bis (e.g. Waxler Theory) ------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') filepowerbis
  WRITE (NDSO,*) 'Power coefficient taken from file:',trim(filepowerbis)
! ------- Name of NetCDF wind file------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') filewind
  WRITE (NDSO,*) 'Wind taken from file:',trim(filewind)
! ------- generic name for summarize at station------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') filestagen
  WRITE (NDSO,*) 'station summation starting with:',trim(filestagen)
! ------- Maximum distance allowed------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,*) MAXDISTANCE
  WRITE (NDSO,*) 'Maximum distance of sources:',MAXDISTANCE 
!
! Input for seismic stations
!
! ------- Nb of stations------------

  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,*) NSTATIONS
  WRITE (NDSO,*) 'Number of seismic stations:',NSTATIONS
  ! Define the size (allocate) of STATIONS, LONO and LATO
! Name of stations need to be written 
  ALLOCATE(STATIONS(NSTATIONS),LONO(NSTATIONS),LATO(NSTATIONS))
  STATIONS(:)=''
  DO I=1,NSTATIONS
    CALL NEXTLN(COMSTR, NDSI, NDSE)
    READ (NDSI,'(A)') STATIONS(I)
  END DO
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,*) LONO(:),LATO(:)
  DO I=1,NSTATIONS
     WRITE (NDSO,*) 'Station:',I,' ',TRIM(STATIONS(I)),', LON = ',LONO(I),', LAT=',LATO(I)
  END DO
! ------- Name of 2nd order pressure variable in WW3 files------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') varp2l
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,*) NTH, TH1
! ------- Path name for Storage of Results ------------
  CALL NEXTLN(COMSTR, NDSI, NDSE)
  READ (NDSI,'(A)') pathNCOUT
  WRITE (NDSO,*) 'Path to storage :',trim(pathNCOUT)

WRITE(6,*) 'COUCOU P2L'
!
! -------------------------------------------
! 1.b.2. Read file with power coefficient
! -------------------------------------------
iret = nf90_open(filepower, NF90_NOWRITE, ncid)
  IF(iret.NE.0) THEN
    WRITE(NDSO,'(2A)') 'ERROR : no file :',trim(filepower)
    STOP
  ELSE
    WRITE(6,'(2A)') 'File opened :',trim(filepower)
  END IF

! Only for regular grid
! 
! ------- lat and lon dimensions------------
iret=nf90_inq_dimid(ncid,'y',dimid(2))
call check_err(iret,'Error in nf90_inq_dimid latitude')
iret=nf90_inquire_dimension(ncid,dimid(2),len=dimln(2))
nlatmod=dimln(2)
iret=nf90_inq_dimid(ncid,'x',dimid(3))
call check_err(iret,'Error in nf90_inq_dimid longitude')
iret=nf90_inquire_dimension(ncid,dimid(3),len=dimln(3))
nlonmod=dimln(3)
iret=nf90_inq_dimid(ncid,'f',dimid(1))
call check_err(iret,'Error in nf90_inq_dimid frequency')
iret=nf90_inquire_dimension(ncid,dimid(1),len=dimln(1))
nfmod=dimln(1)

! ------- model variable latitude------------
IF(ALLOCATED(modlat)) deallocate(modlat)
allocate(modlat(nlatmod))
iret=nf90_inq_varid(ncid,'latitude',varidtmp)
call check_err(iret,'Error in nf90_inq_varid latitude')
iret=nf90_get_var(ncid, varidtmp,modlat(:))
call check_err(iret,'Error in nf90_get_var latitude')
latmin=minval(modlat)
latmax=maxval(modlat)

! ------- model variable longitude------------
IF(ALLOCATED(modlon)) deallocate(modlon)
allocate(modlon(nlonmod))
iret=nf90_inq_varid(ncid,'longitude',varidtmp)
call check_err(iret,'Error in nf90_inq_varid longitude')
iret=nf90_get_var(ncid, varidtmp,modlon(:))
call check_err(iret,'Error in nf90_get_var longitude')
lonmin=minval(modlon)
lonmax=maxval(modlon)

! ------- model variable longitude------------
IF(ALLOCATED(modfreq)) deallocate(modfreq)
allocate(modfreq(nfmod))
iret=nf90_inq_varid(ncid,'frequency',varidtmp)
call check_err(iret,'Error in nf90_inq_varid frequency')
iret=nf90_get_var(ncid, varidtmp,modfreq(:))
call check_err(iret,'Error in nf90_get_var frequency')

! ------- 2D pointer lat and lon------------
IF(ASSOCIATED(alat)) deallocate(alat)
allocate(alat(nlonmod,nlatmod))
IF(ASSOCIATED(alon)) deallocate(alon)
allocate(alon(nlonmod,nlatmod))
DO i=1,nlonmod
  alat(i,:)=modlat(:)
END DO
DO i=1,nlatmod
  alon(:,i)=modlon(:)
END DO

! ------- model variable Power coefficient------------
IF (ALLOCATED(PMAP)) DEALLOCATE(PMAP)
ALLOCATE(PMAP(NLONMOD,NLATMOD,NFMOD))
iret=nf90_inq_varid(ncid,'pmap',varidtmp)
iret=nf90_get_var(ncid,varidtmp,pmap(1:NLONMOD,1:NLATMOD,1:NFMOD))

iret=nf90_close(ncid)

iret = nf90_open(filepowerbis, NF90_NOWRITE, ncid)
  IF(iret.NE.0) THEN
    WRITE(NDSO,'(3A)') 'No file :',trim(filepowerbis),' then nothing is done for Power bis'
    FLAG_POWBIS = .FALSE.
  ELSE
    WRITE(6,'(2A)') 'File opened :',trim(filepowerbis)
    FLAG_POWBIS =.TRUE.
    IF (ALLOCATED(PMAPBIS)) DEALLOCATE(PMAPBIS)
    ALLOCATE(PMAPBIS(NLONMOD,NLATMOD,NFMOD))
    iret=nf90_inq_varid(ncid,'pmap',varidtmp)
    iret=nf90_get_var(ncid,varidtmp,pmapbis(1:NLONMOD,1:NLATMOD,1:NFMOD))

    iret=nf90_close(ncid)
  END IF


! SOIT on rentre le calcule de Pmap ? Soit non (je propose que non) DONC On calcule juste l'atténuation par case et par azimuth
! dpt devient pmap 

nt=floor((REAL(JDAY2-JDAY1)*24+DATE2(4)-DATE1(4))/(dt*24))+1 ! Number of time steps
ALLOCATE(TIMES(nt))



DO I=1,NSTATIONS

! -----------------------------------------------------------------------------
! 1.b.3.Read Wind file for station /!\ IF-LOOP to handle any difference in time steps
! -----------------------------------------------------------------------------
! ------- filewind=fileswind(I)------------

  iret = nf90_open(filewind, NF90_NOWRITE, nciw)
  IF(iret.NE.0) THEN
    WRITE(NDSO,'(2A)') 'ERROR : no file :',trim(filewind)
    STOP
  END IF

! ------- azimuth and time dimensions------------
  iret=nf90_inq_dimid(nciw,'a',dimid(1))
  call check_err(iret,'Error in nf90_inq_dimid azimuth')
  iret=nf90_inquire_dimension(nciw,dimid(1),len=dimln(1))
  nwindazi=dimln(1)
  iret=nf90_inq_dimid(nciw,'t',dimid(2))
  call check_err(iret,'Error in nf90_inq_dimid timewind')
  iret=nf90_inquire_dimension(nciw,dimid(2),len=dimln(2))
  ntimewind=dimln(2)

! ------- wind azimut variable ------------
  IF(ALLOCATED(aziwind)) deallocate(aziwind)
  allocate(aziwind(nwindazi))
  iret=nf90_inq_varid(nciw,'azimuth',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid azimuth')
  iret=nf90_get_var(nciw, varidtmp,aziwind(:))
  call check_err(iret,'Error in nf90_get_var azimuth')

! ------- wind time variable------------
  IF(ALLOCATED(timewind)) deallocate(timewind)
  allocate(timewind(ntimewind))
  iret=nf90_inq_varid(nciw,'time',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid timewind')
  iret=nf90_get_var(nciw, varidtmp,timewind(:))
  call check_err(iret,'Error in nf90_get_var timewind')
  timewind(:) = timewind(:) - 726834
  wtmin=minval(timewind)
  wtmax=maxval(timewind)

  IF(wtmin.GT.JDAY1) THEN
    WRITE(NDSO,'(A)') 'ERROR wind file starts after start date'
    STOP
  END IF
  IF(wtmax.LT.JDAY2) THEN
    WRITE(NDSO,'(A)') 'ERROR wind file ends before end date'
    STOP
  END IF

  dtw=(wtmax-wtmin)/real(ntimewind-1)
!  WRITE(6,*) 'wtmax is :', wtmax
!  WRITE(6,*) 'wtmin is : ', wtmin
  WRITE(6,*) 'ntimewind is:',ntimewind
  IF (DTW.GT.DT) THEN
    DT=DTW
    WRITE(6,*) 'DT has been updated to DTW' ! CHANGE TO BE DONE : option not to update but to duplicate wind
    nt=floor((REAL(JDAY2-JDAY1)*24+DATE2(4)-DATE1(4))/(dt*24))+1 !     Number of time steps
    ALLOCATE(TIMES(nt))
  END IF 
! ------- Vratio variable (azi,time)------------
  IF(ALLOCATED(Vratio)) deallocate(Vratio)
  allocate(Vratio(nwindazi,ntimewind))
  iret=nf90_inq_varid(nciw,'Vratio',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid Vratio')
  iret=nf90_get_var(nciw, varidtmp,Vratio(:,:))
  call check_err(iret,'Error in nf90_get_var timewind')
! .........................................................................
! 
! -----------------------------------------------------------------------------
! 1.b.4.Read station file (for source summation)
! -----------------------------------------------------------------------------
  filesta=trim(filestagen)//trim(STATIONS(I))//'.nc'
  iret = nf90_open(filesta, NF90_NOWRITE, ncis)
  IF(iret.NE.0) THEN
    WRITE(NDSO,'(2A)') 'ERROR : no file :',trim(filesta)
    STOP
  END IF
! ------- Dimensions------------
  iret=nf90_inq_dimid(ncis,'az',dimid(1))
  call check_err(iret,'Error in nf90_inq_dimid azimuth')
  iret=nf90_inquire_dimension(ncis,dimid(1),len=dimln(1))
  nstaazi=dimln(1)
! 
  iret=nf90_inq_dimid(ncis,'dim1',dimid(2))
  call check_err(iret,'Error in nf90_inq_dimid ilon')
  iret=nf90_inquire_dimension(ncis,dimid(2),len=dimln(2))
  nilon=dimln(2)
  nilat=dimln(2)
  nweight=dimln(2)
  ndis=dimln(2)
! 
  iret=nf90_inq_dimid(ncis,'az',dimid(3))
  call check_err(iret,'Error in nf90_inq_dimid startazi')
  iret=nf90_inquire_dimension(ncis,dimid(3),len=dimln(3))
  nstartazi=dimln(3)
  nendazi=dimln(3)
! ------- sta azimut variable------------ 
  IF(ALLOCATED(azista)) deallocate(azista)
  allocate(azista(nstaazi))
  iret=nf90_inq_varid(ncis,'azimuth',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid azimuth')
  iret=nf90_get_var(ncis, varidtmp,azista(:))
  call check_err(iret,'Error in nf90_get_var azimuth')
! ------- sta ilon variable------------ 
  IF(ALLOCATED(ilon)) deallocate(ilon)
  allocate(ilon(nilon))
  iret=nf90_inq_varid(ncis,'ilon',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid ilon')
  iret=nf90_get_var(ncis, varidtmp,ilon(:))
  call check_err(iret,'Error in nf90_get_var ilon')
! ------- sta ilat variable------------ 
  IF(ALLOCATED(ilat)) deallocate(ilat)
  allocate(ilat(nilat))
  iret=nf90_inq_varid(ncis,'ilat',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid ilat')
  iret=nf90_get_var(ncis, varidtmp,ilat(:))
  call check_err(iret,'Error in nf90_get_var ilat')
! ------- sta weight variable------------ 
  IF(ALLOCATED(weight)) deallocate(weight)
  allocate(weight(nweight))
  iret=nf90_inq_varid(ncis,'weight',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid weight')
  iret=nf90_get_var(ncis, varidtmp,weight(:))
  call check_err(iret,'Error in nf90_get_var weight')
! ------- sta distance variable------------ 
  IF(ALLOCATED(DISTANCE)) deallocate(DISTANCE)
  allocate(DISTANCE(ndis))
  iret=nf90_inq_varid(ncis,'distance',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid distance')
  iret=nf90_get_var(ncis, varidtmp,DISTANCE(:))
  call check_err(iret,'Error in nf90_get_var distance')
! ------- sta startazi variable------------ 
  IF(ALLOCATED(startazi)) deallocate(startazi)
  allocate(startazi(nstartazi))
  iret=nf90_inq_varid(ncis,'startazi',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid startazi')
  iret=nf90_get_var(ncis, varidtmp,startazi(:))
  call check_err(iret,'Error in nf90_get_var startazi')
! ------- sta endazi variable------------ 
  IF(ALLOCATED(endazi)) deallocate(endazi)
  allocate(endazi(nendazi))
  iret=nf90_inq_varid(ncis,'endazi',varidtmp)
  call check_err(iret,'Error in nf90_inq_varid endazi')
  iret=nf90_get_var(ncis, varidtmp,endazi(:))
  call check_err(iret,'Error in nf90_get_var endazi')
! IF ANY Ilat >nlatmod... do something (like open lat/lon and see the differences with model? or stop and rerun the model)

  IF (NTH.EQ.0) THEN
    DTH = (360)/(nstaazi-1)
  ELSE
    DTH = 360/NTH
  END IF
WRITE(6,*) 'DTH is ',DTH

! -----------------------------------------------------------------------------
! 1.b.5.Open Source P2L (initialize)
! -----------------------------------------------------------------------------
  WRITE(filenc,'(A,I4,I2.2,A)') trim(pathp2l),DATE1(1:2),'_p2l.nc'
  iret=nf90_open(filenc,NF90_NOWRITE,ncip)
  IF(iret.NE.0) THEN
    WRITE(NDSO,'(2A)') 'ERROR : no file :',trim(filenc)
    STOP
  ELSE
    WRITE(NDSO,'(2A)') 'File opened with nf90_open (ncip) :',trim(filenc)
  END IF
  iret=nf90_inq_dimid(ncip,'time',dimt)
  iret=nf90_inquire_dimension(ncip,dimt,len=NTMOD)
  iret=nf90_inq_dimid(ncip,'f',dimf)
  iret=nf90_inquire_dimension(ncip,dimf,len=NK)
  
  iret=nf90_inq_dimid(ncip,'latitude',dimf)
  iret=nf90_inquire_dimension(ncip,dimf,len=NLATNC)
  iret=nf90_inq_dimid(ncip,'longitude',dimf)
  iret=nf90_inquire_dimension(ncip,dimf,len=NLONNC)
  
  ALLOCATE(TIMESMOD(NTMOD))
  IF (ALLOCATED(FREQ1)) DEALLOCATE(FREQ1)
  IF (ALLOCATED(FREQ2)) DEALLOCATE(FREQ2)
  IF (ALLOCATED(DF)) DEALLOCATE(DF)
  ALLOCATE(FREQ(NK),FREQ1(NK),FREQ2(NK),DF(NK))
  iret=nf90_inq_varid(ncip,'f',varf)
  iret=nf90_get_var(ncip,varf,FREQ)
  iret=nf90_inq_varid(ncip,'time',vartime)
  iret=nf90_get_var(ncip,vartime,TIMESMOD)
  dtmod=(TIMESMOD(NTMOD)-TIMESMOD(1))/real(ntmod-1)
  WRITE(6,*) 'DTMOD is : ', DTMOD
  FREQ=2*FREQ
  FREQ1(1)=FREQ(1)
  DO K=2,NK
    FREQ1(K)=0.5*(FREQ(K)+FREQ(K-1))
    FREQ2(K-1)=FREQ1(K)
  END DO
  FREQ2(NK)=FREQ(NK)
  DF(:)=FREQ2(:)-FREQ1(:)


! -----------------------------------------------------------------------------
! Deal with problem if size WW3.nc not equal to the size of coeff
! -----------------------------------------------------------------------------
! ------- DT issues ------------ 
  IF (DT.NE.DTMOD) THEN
    STEPMOD = floor(DT/DTMOD)
    IF (mod(DT,dtmod).NE.0) THEN
      WRITE(6,*) 'Error DT should be a multiple of dtmod : ',dtmod
      STOP
    END IF
  END IF
! ------- LON issues ------------
  IF (NLONNC.NE.NLONMOD) THEN
    WRITE(6,*) 'Not same size of longitude between model NC and Power coeff'
    STOP    
  END IF
! ------- LAT issues ------------
! Write(6,*) NLATNC, NLATMOD
  IF (NLATNC.NE.NLATMOD) THEN
    IF (NLATNC.GT.NLATMOD) THEN
      NLATNC=NLATMOD
    ELSE
      WRITE(6,*) 'Not same size of latitude between model NC and Power coeff'
      STOP 
    END IF  
  END IF
! ------- Frequency issues ------------
  IF (NK.NE.NFMOD) THEN
    WRITE(6,*) 'Warning: Not same size of frequency between model NC and Power coeff'
    
    IF ((NK.GT.NFMOD).AND.(ANY(FREQ.EQ.MODFREQ(1)))) THEN
      indfreqdeb = minloc(FREQ,dim=1, mask=FREQ.EQ.MODFREQ(1))
      ! WRITE(6,*) 'Ind =', SIZE(FREQ(indfreqdeb:indfreqdeb+NFMOD-1))
      IF (ALLOCATED(FREQB)) DEALLOCATE(FREQB)
      IF (NK.GE.indfreqdeb+NFMOD-1) THEN
        ALLOCATE(FREQB(NFMOD))!-1 added
        FREQB = FREQ(indfreqdeb:indfreqdeb+NFMOD-1)!-1 added
        IF (ALL(FREQB-MODFREQ.EQ.0)) THEN
          IF (ALLOCATED(FREQ)) DEALLOCATE(FREQ)
          ALLOCATE(FREQ(NFMOD))
          FREQ = FREQB
          FREQB = FREQ1(indfreqdeb:indfreqdeb+NFMOD-1)
          IF (ALLOCATED(FREQ1)) DEALLOCATE(FREQ1)
          ALLOCATE(FREQ1(NFMOD))
          FREQ1 = FREQB
          FREQB = FREQ2(indfreqdeb:indfreqdeb+NFMOD-1)
          IF (ALLOCATED(FREQ2)) DEALLOCATE(FREQ2)
          ALLOCATE(FREQ2(NFMOD))
          FREQ2 = FREQB
          NK = NFMOD
          WRITE(6,*) 'Problem of frequency solved (frequencies p2l contains frequencies of Power)'
        ELSE
          WRITE(6,*) 'Error : Frequencies do not match'
          STOP
        END IF
      ELSE
        IF (ALLOCATED(MODFREQB)) DEALLOCATE(MODFREQB)
        ALLOCATE(MODFREQB(NK+1-indfreqdeb))
        MODFREQB = MODFREQ(1:NK-1+indfreqdeb)
        ALLOCATE(FREQB(NK+1-indfreqdeb))
        FREQB = FREQ(indfreqdeb:NK)
        IF (ALL(FREQB.EQ.MODFREQB)) THEN
          FREQ = FREQB
          FREQB = FREQ1(indfreqdeb:NK)
          IF (ALLOCATED(FREQ1)) DEALLOCATE(FREQ1)
          ALLOCATE(FREQ1(NK+1-indfreqdeb))
          FREQ1 = FREQB
          FREQB = FREQ2(indfreqdeb:NK)
          IF (ALLOCATED(FREQ2)) DEALLOCATE(FREQ2)
          ALLOCATE(FREQ2(NK+1-indfreqdeb))
          FREQ2 = FREQB
          NK = NK+1-indfreqdeb
          WRITE(6,*) 'Problem of frequency solved (frequencies p2l and frequencies of Power overlap)'
        ELSE
          WRITE(6,*) 'Error : Frequencies do not match'
          STOP
        END IF
      END IF
    END IF
  END IF

  iret=nf90_inq_varid(ncip,varp2l,varidp2l)
  IF(iret.NE.0)THEN
    WRITE(6,*) 'NetCDF variable not found:',TRIM(varp2l),' in file ',TRIM(filenc)
    STOP
  END IF

  iret=nf90_get_att(ncip,varidp2l,'scale_factor',SCALEP2L)
  iret=nf90_get_att(ncip,varidp2l,'units',UNITP2L)
  iret=nf90_get_att(ncip,varidp2l,'_FillValue',FILLP2L)
! -----------------------------------------------------------------------------
! Transform integer array into P2L via inverse'log'
! -----------------------------------------------------------------------------
  WRITE(6,*) SCALEP2L
  imin=-2**15+1
  imax=2**15
  ! NTAB=imax-imin+1
  ALLOCATE(TABP2F(IMIN:IMAX))
  TABP2F=(/(I,I=IMIN,IMAX)/)
  lenu=LEN_TRIM(UNITP2L)
  IF (INDEX(UNITP2L,'log10(m4s+1E-12').EQ.1) &
   tabP2F=exp(lg10*tabP2F*DBLE(SCALEP2L))-(DBLE(1E-12)-DBLE(1E-16))
  IF (INDEX(UNITP2L,'log10(m4s+0.01)').EQ.1) &
   tabP2F=exp(lg10*tabP2F*SCALEP2L)-DBLE(0.009999)
  IF (INDEX(UNITP2L,'log10(Pa2s+1E-12)').EQ.1) &
   tabP2F=exp(lg10*tabP2F*SCALEP2L)-(1E-12-1E-16)
  IF (INDEX(UNITP2L,'log10(Pa s+1E-12)').EQ.1) &
   tabP2F=exp(lg10*tabP2F*SCALEP2L)-(1E-12-1E-16)
  IF (INDEX(UNITP2L,'log10(Pa2 m2 s+1E-12').EQ.1) &
   tabP2F=exp(lg10*tabP2F*SCALEP2L)-((1E-12)-(1E-16))

  tabP2f(imin)=0. ! Should correspond to Fill_value
 ! WRITE(6,*) 'TABP2F is negative for : ',COUNT(TABP2F.LT.0),' values'

  ntrmax=8
  itr=0
  ALLOCATE(P2ft(NLONMOD,NLATMOD,NK,NTRMAX))!,P2FT0(NLONMOD,NLATMOD,NK,NTRMAX))!,P2F(NLONMOD,NLATMOD,NK))
  P2ft(:,:,:,:)=imin
 ! P2ft0(:,:,:,:)=imin
!  IF(ALLOCATED(DISTM1)) DEALLOCATE(DISTM1)
!  IF(ALLOCATED(DISTM2)) DEALLOCATE(DISTM2)
!  IF(ALLOCATED(AMPM1)) DEALLOCATE(AMPM1)
!  IF(ALLOCATED(AMPM2)) DEALLOCATE(AMPM2)
!  IF(ALLOCATED(F_spec_tot1)) DEALLOCATE(F_spec_tot1)
!  IF(ALLOCATED(F_spec_tot2)) DEALLOCATE(F_spec_tot2)
!  IF(ALLOCATED(DISTMBIS1)) DEALLOCATE(DISTMBIS1)
!  IF(ALLOCATED(DISTMBIS2)) DEALLOCATE(DISTMBIS2)
!  IF(ALLOCATED(AMPMBIS1)) DEALLOCATE(AMPMBIS1)
!  IF(ALLOCATED(AMPMBIS2)) DEALLOCATE(AMPMBIS2)
!  IF(ALLOCATED(F_spec_totBIS1)) DEALLOCATE(F_spec_totBIS1)
!  IF(ALLOCATED(F_spec_totBIS2)) DEALLOCATE(F_spec_totBIS2)

  IF(ALLOCATED(F_spec1)) DEALLOCATE(F_spec1)
  IF(ALLOCATED(F_spec2)) DEALLOCATE(F_spec2)

  IF(ALLOCATED(F_specBIS1)) DEALLOCATE(F_specBIS1)
  IF(ALLOCATED(F_specBIS2)) DEALLOCATE(F_specBIS2)

  IF(ALLOCATED(Amp_t_az_Delf1)) DEALLOCATE(Amp_t_az_Delf1)!(:,:,:),
  IF(ALLOCATED(Amp_t_az_Delf2)) DEALLOCATE(Amp_t_az_Delf2)!(:,:,:),
  IF(ALLOCATED(Amp_t_az_DelfBIS1)) DEALLOCATE(Amp_t_az_DelfBIS1)!(:,:,:),
  IF(ALLOCATED(Amp_t_az_DelfBIS2)) DEALLOCATE(Amp_t_az_DelfBIS2)!(:,:,:)
  IF(ALLOCATED(f_t_az_mamp1)) DEALLOCATE(f_t_az_mamp1)
  IF(ALLOCATED(f_t_az_mamp2)) DEALLOCATE(f_t_az_mamp2)
  IF(ALLOCATED(f_t_az_mampBIS1)) DEALLOCATE(f_t_az_mampBIS1)
  IF(ALLOCATED(f_t_az_mampBIS2)) DEALLOCATE(f_t_az_mampBIS2)
! Double precision, allocatable 
  IF(ALLOCATED(Az_t_f_mamp1)) DEALLOCATE(Az_t_f_mamp1)
  IF(ALLOCATED(Az_t_f_mamp2)) DEALLOCATE(Az_t_f_mamp2)
  IF(ALLOCATED(Az_t_f_mampBIS1)) DEALLOCATE(Az_t_f_mampBIS1)
  IF(ALLOCATED(Az_t_f_mampBIS2)) DEALLOCATE(Az_t_f_mampBIS2)
! Double precision, allocatable 
  IF(ALLOCATED(Amp_t_az_sum1)) DEALLOCATE(Amp_t_az_sum1) 
  IF(ALLOCATED(Amp_t_az_sum2)) DEALLOCATE(Amp_t_az_sum2)
  IF(ALLOCATED(Amp_t_az_sumBIS1)) DEALLOCATE(Amp_t_az_sumBIS1)
  IF(ALLOCATED(Amp_t_az_sumBIS2)) DEALLOCATE(Amp_t_az_sumBIS2)
! Double precision, allocatable 
  IF(ALLOCATED(Amp_t_az_mf1)) DEALLOCATE(Amp_t_az_mf1) 
  IF(ALLOCATED(Amp_t_az_mf2)) DEALLOCATE(Amp_t_az_mf2)
  IF(ALLOCATED(Amp_t_az_mfBIS1)) DEALLOCATE(Amp_t_az_mfBIS1)
  IF(ALLOCATED(Amp_t_az_mfBIS2)) DEALLOCATE(Amp_t_az_mfBIS2)
! Double precision, allocatable 
  IF(ALLOCATED(Amp_t_f_sum1)) DEALLOCATE(Amp_t_f_sum1)
  IF(ALLOCATED(Amp_t_f_sum2))   DEALLOCATE(Amp_t_f_sum2)
  IF(ALLOCATED(Amp_t_f_sumBIS1)) DEALLOCATE(Amp_t_f_sumBIS1)
  IF(ALLOCATED(Amp_t_f_sumBIS2)) DEALLOCATE(Amp_t_f_sumBIS2)
! Double precision, allocatable 
  IF(ALLOCATED(Amp_t_f_maz1)) DEALLOCATE(Amp_t_f_maz1) 
  IF(ALLOCATED(Amp_t_f_maz2)) DEALLOCATE(Amp_t_f_maz2)
  IF(ALLOCATED(Amp_t_f_mazBIS1)) DEALLOCATE(Amp_t_f_mazBIS1)
  IF(ALLOCATED(Amp_t_f_mazBIS2)) DEALLOCATE(Amp_t_f_mazBIS2)

  ALLOCATE(F_spec1(NSTAAZI,NK,NT),Amp_t_az_Delf1(NSTAAZI,5,NT))
  ALLOCATE(F_spec2(NSTAAZI,NK,NT),Amp_t_az_Delf2(NSTAAZI,5,NT))
  ALLOCATE(Amp_t_f_sum1(NK,NT),Az_t_f_mamp1(NK,NT),Amp_t_f_maz1(NK,NT))
  ALLOCATE(Amp_t_f_sum2(NK,NT),Az_t_f_mamp2(NK,NT),Amp_t_f_maz2(NK,NT))
  ALLOCATE(Amp_t_az_sum1(NK,NT),f_t_az_mamp1(NK,NT),Amp_t_az_mf1(NK,NT))
  ALLOCATE(Amp_t_az_sum2(NK,NT),f_t_az_mamp2(NK,NT),Amp_t_az_mf2(NK,NT))

!  ALLOCATE(DISTM1(NSTAAZI,NK,NT),AMPM1(NSTAAZI,NK,NT),F_spec1(NSTAAZI,NK,NT))
!  ALLOCATE(DISTM2(NSTAAZI,NK,NT),AMPM2(NSTAAZI,NK,NT),F_spec2(NSTAAZI,NK,NT))
!  ALLOCATE(F_spec_tot1(NK,NT))
!  ALLOCATE(F_spec_tot2(NK,NT))  

  IF(FLAG_POWBIS)THEN
!    ALLOCATE(DISTMBIS1(NSTAAZI,NK,NT),AMPMBIS1(NSTAAZI,NK,NT),F_specBIS1(NSTAAZI,NK,NT))
!    ALLOCATE(F_spec_totBIS1(NK,NT)) 
!    ALLOCATE(DISTMBIS2(NSTAAZI,NK,NT),AMPMBIS2(NSTAAZI,NK,NT),F_specBIS2(NSTAAZI,NK,NT))
!    ALLOCATE(F_spec_totBIS2(NK,NT)) 
    ALLOCATE(F_specBIS1(NSTAAZI,NK,NT),Amp_t_az_DelfBIS1(NSTAAZI,5,NT))
    ALLOCATE(F_specBIS2(NSTAAZI,NK,NT),Amp_t_az_DelfBIS2(NSTAAZI,5,NT))
    ALLOCATE(Amp_t_f_sumBIS1(NK,NT),Az_t_f_mampBIS1(NK,NT),Amp_t_f_mazBIS1(NK,NT))
    ALLOCATE(Amp_t_f_sumBIS2(NK,NT),Az_t_f_mampBIS2(NK,NT),Amp_t_f_mazBIS2(NK,NT))
    ALLOCATE(Amp_t_az_sumBIS1(NK,NT),f_t_az_mampBIS1(NK,NT),Amp_t_az_mfBIS1(NK,NT))
    ALLOCATE(Amp_t_az_sumBIS2(NK,NT),f_t_az_mampBIS2(NK,NT),Amp_t_az_mfBIS2(NK,NT))
  END IF

!  IF (ALLOCATED(FREQ1)) DEALLOCATE(FREQ1)
!  IF (ALLOCATED(FREQ2)) DEALLOCATE(FREQ2)
!  IF (ALLOCATED(DF)) DEALLOCATE(DF)
!  ALLOCATE(FREQ1(NK),FREQ2(NK),DF(NK))
!  FREQ1(1)=FREQ(1)
!  DO K=2,NK
!    FREQ1(K)=0.5*(FREQ(K)+FREQ(K-1))
!    FREQ2(K-1)=FREQ1(K)
!  END DO
!  FREQ2(NK)=FREQ(NK)
!  DF(:)=FREQ2(:)-FREQ1(:)

! -----------------------------------------------------------------------------
! 1.b.6. Initialize NC file out (1 per station)
! -----------------------------------------------------------------------------

  dimlno(1) = NT
  dimlno(2) = 1
  dimlno(3) = 16
  dimlno(4) = NK
  dimlno(5) = NSTAAZI
  dimlno(6) = 5
IF (MAXDISTANCE.GE.10000)THEN
  WRITE(filencout1,'(A,A,A,I4,I2.2,A,I4,I2.2,A,A,A,I5,A)') trim(pathNCOUT),trim(stations(I)),'_',DATE1(1),DATE1(2),&
'-',DATE2(1),DATE2(2),'_',trim(REFTAG),'_WindSta_distmax',int(MAXDISTANCE),'.nc'
  WRITE(filencout2,'(A,A,A,I4,I2.2,A,I4,I2.2,A,A,A,I5,A)') trim(pathNCOUT),trim(stations(I)),'_',DATE1(1),DATE1(2),&
'-',DATE2(1),DATE2(2),'_',trim(REFTAG),'_NoWind_distmax',int(MAXDISTANCE),'.nc'
ELSE
  WRITE(filencout1,'(A,A,A,I4,I2.2,A,I4,I2.2,A,A,A,I4,A)') trim(pathNCOUT),trim(stations(I)),'_',DATE1(1),DATE1(2),&
'-',DATE2(1),DATE2(2),'_',trim(REFTAG),'_WindSta_distmax',int(MAXDISTANCE),'.nc'
  WRITE(filencout2,'(A,A,A,I4,I2.2,A,I4,I2.2,A,A,A,I4,A)') trim(pathNCOUT),trim(stations(I)),'_',DATE1(1),DATE1(2),&
'-',DATE2(1),DATE2(2),'_',trim(REFTAG),'_NoWind_distmax',int(MAXDISTANCE),'.nc'
END IF


!  WRITE(filencout,'(A,A,I4,I2.2,A,I4,I2.2,A,A,A,A,A)')  trim(stations(I)),'_',&
!     DATE1(1),DATE1(2),'-',DATE2(1),DATE2(2),'_',trim(REFTAG),'_',trim(WINDTAG),'.nc'
  WRITE(6,*) 'Initializing output file:',trim(filencout1)
  CALL W3CRNC (4,filencout1,ncio1,dimido,dimlno,varido1,FLAG_POWBIS)
  !iret = nf90_put_var(ncio1,varido(27),I,1)
  iret = nf90_put_var(ncio1,varido1(2),trim(stations(I)),start=(/1,1/),count=(/LEN_TRIM(STATIONS(I)),1/))
  iret = nf90_put_var(ncio1,varido1(6),FREQ(1:NK)/2)
  iret = nf90_put_var(ncio1,varido1(7),FREQ1(1:NK))
  iret = nf90_put_var(ncio1,varido1(8),FREQ2(1:NK))
  iret = nf90_put_var(ncio1,varido1(9),AZISTA(1:NSTAAZI))
  iret = nf90_put_var(ncio1,varido1(26),(/0.1,0.2,0.3,0.4,0.5/))
  WRITE(6,*) 'Initializing output file:',trim(filencout2)
  CALL W3CRNC (4,filencout2,ncio2,dimido,dimlno,varido2,FLAG_POWBIS)
  !iret = nf90_put_var(ncio2,varido(27),I,1)
  iret = nf90_put_var(ncio2,varido2(2),trim(stations(I)),start=(/1,1/),count=(/LEN_TRIM(STATIONS(I)),1/))
  iret = nf90_put_var(ncio2,varido2(6),FREQ(1:NK)/2)
  iret = nf90_put_var(ncio2,varido2(7),FREQ1(1:NK))
  iret = nf90_put_var(ncio2,varido2(8),FREQ2(1:NK))
  iret = nf90_put_var(ncio2,varido2(9),AZISTA(1:NSTAAZI))
  iret = nf90_put_var(ncio2,varido2(26),(/0.1,0.2,0.3,0.4,0.5/))
!===============================================================================  
! 2. TIME-LOOP STARTING
!===============================================================================  
  DO it=1,NT 
    !===============================================================================  
    ! 2.a Open New WW3 model file if needed and get 8 time steps
    !===============================================================================   
    NOW=(IT-1)*DT
    CALL CALDAT(JDAY1+JDAYREF+FLOOR(NOW),id,mm,iyyy)
    WRITE(6,*) 'id/mm/iyyy = ', id,'/',mm,'/',iyyy
    FILENCOLD=FILENC
    WRITE(filenc,'(A,I4,I2.2,A)') TRIM(pathp2l),iyyy,mm,'_p2l.nc'
    WRITE(6,*) 'try to open :', trim(filenc)
    IF (TRIM(FILENC).NE.TRIM(FILENCOLD)) THEN
      iret=nf90_close(ncip)
      DEALLOCATE(TIMESMOD)
      WRITE(6,*) 'New model file:',TRIM(filenc)
      iret=nf90_open(filenc,NF90_NOWRITE,ncip)
      iret=nf90_inq_dimid(ncip,'time',dimt)
      iret=nf90_inquire_dimension(ncip,dimt,len=NTMOD)
      ALLOCATE(TIMESMOD(NTMOD))
      iret=nf90_inq_varid(ncip,'time',vartime)
      iret=nf90_get_var(ncip,vartime,TIMESMOD)
      dtmod=(TIMESMOD(NTMOD)-TIMESMOD(1))/real(ntmod-1)
      iret=nf90_inq_varid(ncip,varp2l,varidp2l)
      iret=nf90_get_att(ncip,varidp2l,'scale_factor',SCALEP2LTMP)
      iret=nf90_get_att(ncip,varidp2l,'units',UNITP2LTMP)
      iret=nf90_get_att(ncip,varidp2l,'_FillValue',FILLP2LTMP)
      IF (UNITP2LTMP(1:min(len_trim(UNITP2L),len_trim(UNITP2LTMP))).NE.&
UNITP2L(1:min(len_trim(UNITP2L),len_trim(UNITP2LTMP)))) THEN
        WRITE(6,*) 'ERROR: CHANGE IN UNITS IN NEW *.nc FILE'
        STOP
      END IF
      IF (SCALEP2LTMP.NE.SCALEP2L) THEN
        WRITE(6,*) 'ERROR: CHANGE IN SCALE FACTOR IN NEW *.nc FILE'
        STOP
      END IF
      IF (FILLP2LTMP.NE.FILLP2L) THEN 
        WRITE(6,*) 'ERROR: CHANGE IN FILL VAL IN NEW *.nc FILE'
        STOP
      END IF

    END IF

    IF ((mod(it,ntrmax).EQ.1)) THEN
      nti=1+(TIMESMOD(NTMOD)-JDAY1)/DT ! HERE NTI is the limit of the source file .nc
      ntr=min(ntrmax,min(nt,nti)-it+1)
      WRITE(6,*) 'NTI is ', NTI,'-------- NTR is :', NTR
      ITIME=1+FLOOR((JDAY1+NOW-TIMESMOD(1))/DTMOD)
      WRITE(6,*) 'NTI is ', NTI
      ! ------- Read for each calculated time steps ------------
      IF (DTMOD.NE.DT) THEN
        WRITE(6,*) 'Before nf90 get var (If DTMOD.NE.DT)'
        WRITE(6,*) 'Size of (lon,lat, freq,Time) requested',NLONMOD,NLATMOD,NK,NTR
        WRITE(6,*) 'Start time index and step time:',ITIME,STEPMOD
        iret=nf90_get_var(ncip,varidp2l,P2FT(1:NLONMOD,1:NLATMOD,1:NK,1:NTR), &
                      start=(/1,1,1,ITIME/),count=(/NLONMOD,NLATMOD,NK,NTR/))
        WRITE(6,*) 'After nf90 get var (IF DTMOD.NE.DT)'
        CALL check_err(iret, 'Error in nf90_get_var varidP2l')
     !   WRITE(6,*) 'If 1 before nf90 inq var time'
        iret=nf90_inq_varid(ncip,'time',vartime)
     !   WRITE(6,*) 'If 1 after nf90 inq var time'
      ELSE
     !   WRITE(6,*) 'If 2 before nf90 get var'
        iret=nf90_get_var(ncip,varidp2l,P2FT(1:NLONMOD,1:NLATMOD,1:NK,1:NTR), &
                      start=(/1,1,1,ITIME/),count=(/NLONMOD,NLATMOD,NK,NTR/))
     !   WRITE(6,*) 'If 2 after nf90 get var'
        CALL check_err(iret, 'Error in nf90_get_var varidP2l')
      END IF
      WRITE(6,'(A,I4,A,I4,A,I2.2,A,I2.2,A,A)') 'Reading NetCDF time step ',IT,' out of ',NT,'  for date ',ID,'/',mm, &
                          ' in file ', trim(filenc)
      IF ((NTR.LT.NTRMAX).AND.(NT.GT.NTI)) THEN
        FILENCOLD = FILENC
     ! ------- READ next source ... ------------
        NOWNEW=(NTI+1-1)*DT
        CALL CALDAT(JDAY1+JDAYREF+FLOOR(NOWNEW),id,mm,iyyy)
        WRITE(filenc,'(A,I4,I2.2,A)') TRIM(pathp2l),iyyy,mm,'_p2l.nc'
        iret=nf90_close(ncip)
        DEALLOCATE(TIMESMOD)
        WRITE(6,*) 'New model file:',TRIM(filenc)
        iret=nf90_open(filenc,NF90_NOWRITE,ncip)
        iret=nf90_inq_dimid(ncip,'time',dimt)
        iret=nf90_inquire_dimension(ncip,dimt,len=NTMOD)
        ALLOCATE(TIMESMOD(NTMOD))
        iret=nf90_inq_varid(ncip,'time',vartime)
        iret=nf90_get_var(ncip,vartime,TIMESMOD)
        dtmod=(TIMESMOD(NTMOD)-TIMESMOD(1))/real(ntmod-1)
        iret=nf90_inq_varid(ncip,varp2l,varidp2l)
        iret=nf90_get_att(ncip,varidp2l,'scale_factor',SCALEP2LTMP)
        iret=nf90_get_att(ncip,varidp2l,'units',UNITP2LTMP)
        iret=nf90_get_att(ncip,varidp2l,'_FillValue',FILLP2LTMP)
        IF (UNITP2LTMP(1:min(len_trim(UNITP2L),len_trim(UNITP2LTMP)))&
.NE.UNITP2L(1:min(len_trim(UNITP2L),len_trim(UNITP2LTMP)))) THEN
        WRITE(6,*) 'ERROR: CHANGE IN UNITS IN NEW *.nc FILE'
        STOP
      END IF
      IF (SCALEP2LTMP.NE.SCALEP2L) THEN
        WRITE(6,*) 'ERROR: CHANGE IN SCALE FACTOR IN NEW *.nc FILE'
        STOP
      END IF
      IF (FILLP2LTMP.NE.FILLP2L) THEN 
        WRITE(6,*) 'ERROR: CHANGE IN FILL VAL IN NEW *.nc FILE'
        STOP
      END IF
        nti=min(ntrmax,nt-it+1)
        ITIME=1+FLOOR((JDAY1+NOWNEW-TIMESMOD(1))/DTMOD)
        ntr = nti-ITIME +1
      !  ------- Read for each calculated time steps ------------
        IF (DTMOD.NE.DT) THEN
          iret=nf90_get_var(ncip,varidp2l,P2FT(1:NLONMOD,1:NLATMOD,1:NK,1:NTR), &
                      start=(/1,1,1,ITIME/),count=(/NLONMOD,NLATMOD,NK,NTR/))
          CALL check_err(iret, 'Error in nf90_get_var varidP2l')
          iret=nf90_inq_varid(ncip,'time',vartime)
        ELSE
          iret=nf90_get_var(ncip,varidp2l,P2FT(1:NLONMOD,1:NLATMOD,1:NK,1:NTR), &
                        start=(/1,1,1,ITIME/),count=(/NLONMOD,NLATMOD,NK,NTR/))
          CALL check_err(iret, 'Error in nf90_get_var varidP2l')
        END IF
        FILENC = FILENCOLD
      END IF ! END of : 'complete with next source'
      itr=1
   !   WRITE(6,*) P2FT(NLONMOD,NLATMOD,1,1), P2FT(NLONMOD,NLATMOD,NK,NTR), FILLP2L
    ELSE
      itr=itr+1
    END IF
    times(it)=now+DBLE(JDAY1);
    !===============================================================================  
    ! 2.b Frequency-loop and Azimuth Loop
    ! azista(nstaazi) (1:1:180 par ex) ilon(nilon)ilat(nilat)weight(nweight) (same length) startazi end azi (length nstaazi)
    !   => two-by-two ways for calculating Vratio : from source or station , with AZE or AZS
    !===============================================================================  
    DO K=1,NK
      DO AZ=1,NSTAAZI
        azi=azista(AZ)
        STARTA=startazi(AZ)
        ENDSTA=endazi(AZ)
!  ------- Vratio at Station : AZE ------------
        Vraze=Vratio(AZ,IT)

	IF (ALLOCATED(DIS)) DEALLOCATE(DIS)
        ALLOCATE(DIS(ENDSTA-STARTA+1))
        DIS=DISTANCE(STARTA:ENDSTA)
        LL=COUNT(DIS.LT.MAXDISTANCE)

        IF (ALLOCATED(WEIGHT2)) DEALLOCATE(WEIGHT2)
        ALLOCATE(WEIGHT2(LL))
        WEIGHT2=PACK(WEIGHT(STARTA:ENDSTA),mask=(DIS.LT.MAXDISTANCE))
        IF (ALLOCATED(ILON0)) DEALLOCATE(ILON0)
        ALLOCATE(ILON0(LL))
        ILON0=PACK(ILON(STARTA:ENDSTA),mask=(DIS.LT.MAXDISTANCE))
        IF (ALLOCATED(ILAT0)) DEALLOCATE(ILAT0)
        ALLOCATE(ILAT0(LL))
        ILAT0=PACK(ILAT(STARTA:ENDSTA),mask=(DIS.LT.MAXDISTANCE))
	IF (ALLOCATED(DIS_MASKED)) DEALLOCATE(DIS_MASKED)
        ALLOCATE(DIS_MASKED(LL))
        DIS_MASKED=PACK(DIS,mask=(DIS.LT.MAXDISTANCE))

        IF (ALLOCATED(SOURCE)) DEALLOCATE(SOURCE)
        IF (ALLOCATED(SOURCEBIS)) DEALLOCATE(SOURCEBIS)
        ALLOCATE(SOURCE(LL))
        IF (ALLOCATED(SOURCEOK)) DEALLOCATE(SOURCEOK)
        ALLOCATE(SOURCEOK(LL))
        IF (ALLOCATED(AMP1)) DEALLOCATE(AMP1)
        ALLOCATE(AMP1(LL))
        IF (ALLOCATED(AMP2)) DEALLOCATE(AMP2)
        ALLOCATE(AMP2(LL))
        IF (ALLOCATED(ATTTMP1)) DEALLOCATE(ATTTMP1)
        ALLOCATE(ATTTMP1(LL))
        IF (ALLOCATED(ATTTMP2)) DEALLOCATE(ATTTMP2)
        ALLOCATE(ATTTMP2(LL))
!  ------- Vratio at Station : AZS (not recommended) ------------
!        IF (ALLOCATED(VRATIOAZ)) DEALLOCATE(VRATIOAZ)
!        ALLOCATE(VRATIOAZ(LL))

        IF(FLAG_POWBIS)THEN
          ALLOCATE(SOURCEBIS(LL))
        END IF
        

        DO PTDIR=1,LL
          IO=ILON0(PTDIR)
          IA=ILAT0(PTDIR)
          SOURCE(PTDIR)=TABP2F(NINT(P2FT(IO,IA,K,min(ITR,NTR))))*PMAP(IO,IA,K)
!  ------- Vratio at Station : AZS (not recommended) ------------
!          CALL DISTAZ(lato(I),lono(I),modlat(ia),modlon(io),aze,azs)
!          AZS2=floor(azs-180)
!          IF (azs2.LE.0) azs2=azs2+360
!          Vratioaz(PTDIR)=Vratio(AZS2,IT) ! This is for Vratio at the station
          IF(FLAG_POWBIS)THEN
            SOURCEBIS(PTDIR)=TABP2F(NINT(P2FT(IO,IA,K,min(ITR,NTR))))*PMAPBIS(IO,IA,K)
          END IF
          IF ((PMAP(IO,IA,K).EQ.0).OR.(P2FT(IO,IA,K,min(ITR,NTR)).EQ.FILLP2L)) THEN
          !  WRITE(6,*) 'HERE'
            SOURCEOK(PTDIR) = .FALSE.
          ELSE
            SOURCEOK(PTDIR) =.TRUE.
          END IF
        END DO ! END DO PTDIR ===========  
!  ------- Vratio at Station : AZS (not recommended) ------------
!        CALL attenuation(DIS_MASKED,FREQ(K),Vratioaz,ATTTMP1,LL) ! Define ATTTMP
!  ------- Vratio at Station : AZE ------------
        CALL attenuationAZE(DIS_MASKED,FREQ(K),Vraze,ATTTMP1,LL) ! Define ATTTMP
        AMP1=ATTTMP1*ATTTMP1*SOURCE
      !  maxiamp=minloc(ATTTMP,dim=1)
!        maxiamp = maxloc(AMP1,dim=1)
!        DISTM1(AZ,K,it)=DIS_MASKED(maxiamp)
!        AMPM1(AZ,K,it)=AMP1(maxiamp)
        F_spec1(AZ,K,it)=SUM(AMP1*WEIGHT2,mask=SOURCEOK)
! ---------------- For Propa without wind
        CALL attenuationNoWind(DIS_MASKED,FREQ(K),ATTTMP2,LL) ! Define ATTTMP
        AMP2=ATTTMP2*ATTTMP2*SOURCE
      !  maxiamp=minloc(ATTTMP,dim=1)
!        maxiamp = maxloc(AMP2,dim=1)
!        DISTM2(AZ,K,it)=DIS_MASKED(maxiamp)
!        AMPM2(AZ,K,it)=AMP2(maxiamp)
        F_spec2(AZ,K,it)=SUM(AMP2*WEIGHT2,mask=SOURCEOK)

        IF(ANY((SOURCE.LT.0).AND.(SOURCEOK)))THEN
         WRITE(6,*) 'THERE are ', COUNT((SOURCE.LT.0).AND.(SOURCEOK)),' negative sources for az = ', AZ,' and k = ', K
        END IF
        IF(FLAG_POWBIS)THEN
          AMP1=ATTTMP1*ATTTMP1*SOURCEBIS
      !  maxiamp=minloc(ATTTMP,dim=1)
!          maxiamp = maxloc(AMP1,dim=1)!*WEIGHT2
!          DISTMBIS1(AZ,K,it)=DIS_MASKED(maxiamp)
!          AMPMBIS1(AZ,K,it)=AMP1(maxiamp)
          F_specBIS1(AZ,K,it)=SUM(AMP1*WEIGHT2,mask=SOURCEOK)
! -------------- For Propa without wind
          AMP2=ATTTMP2*ATTTMP2*SOURCEBIS
!          maxiamp = maxloc(AMP2,dim=1)!*WEIGHT2
!          DISTMBIS2(AZ,K,it)=DIS_MASKED(maxiamp)
!          AMPMBIS2(AZ,K,it)=AMP2(maxiamp)
          F_specBIS2(AZ,K,it)=SUM(AMP2*WEIGHT2,mask=SOURCEOK)
        END IF

       ! WRITE(6,*) DIS(maxiamp)
      END DO !DO AZ=1,NSTAAZI
!      F_spec_tot1(K,it)=SUM(F_spec1(:,K,it)*dth,dim=1)
!      F_spec_tot2(K,it)=SUM(F_spec2(:,K,it)*dth,dim=1)
      Amp_t_f_sum1(K,it)=SUM(F_spec1(:,K,it)*dth,dim=1, mask=(.NOT.ISNAN (F_spec1(:,K,it)*dth))) ! F_spec_tot1
      Amp_t_f_sum2(K,it)=SUM(F_spec2(:,K,it)*dth,dim=1, mask=(.NOT.ISNAN (F_spec2(:,K,it)*dth))) ! F_spec_tot2
      Amp_t_f_maz1(K,it)=maxval(F_spec1(:,K,it),dim=1, mask=(.NOT.ISNAN (F_spec1(:,K,it))))
      Amp_t_f_maz2(K,it)=maxval(F_spec2(:,K,it),dim=1, mask=(.NOT.ISNAN (F_spec2(:,K,it))))
      Az_t_f_mamp1(K,it)=AZISTA(maxloc(F_spec1(:,K,it),dim=1, mask=(.NOT.ISNAN (F_spec1(:,K,it)))))
      Az_t_f_mamp2(K,it)=AZISTA(maxloc(F_spec2(:,K,it),dim=1, mask=(.NOT.ISNAN (F_spec2(:,K,it)))))

      IF(FLAG_POWBIS)THEN
!        F_spec_totBIS1(K,it)=SUM(F_specBIS1(:,K,it)*dth,dim=1)
!        F_spec_totBIS2(K,it)=SUM(F_specBIS2(:,K,it)*dth,dim=1)
        Amp_t_f_sumBIS1(K,it)=SUM(F_specBIS1(:,K,it)*dth,dim=1, mask=(.NOT.ISNAN(F_specBIS1(:,K,it)*dth))) ! F_spec_totBIS1
        Amp_t_f_sumBIS2(K,it)=SUM(F_specBIS2(:,K,it)*dth,dim=1, mask=(.NOT.ISNAN(F_specBIS2(:,K,it)*dth))) ! F_spec_totBIS2
        Amp_t_f_mazBIS1(K,it)=maxval(F_specBIS1(:,K,it),dim=1, mask=(.NOT.ISNAN(F_specBIS1(:,K,it))))
        Amp_t_f_mazBIS2(K,it)=maxval(F_specBIS2(:,K,it),dim=1, mask=(.NOT.ISNAN(F_specBIS2(:,K,it))))
        Az_t_f_mampBIS1(K,it)=AZISTA(maxloc(F_specBIS1(:,K,it),dim=1, mask=(.NOT.ISNAN(F_specBIS1(:,K,it)))))
        Az_t_f_mampBIS2(K,it)=AZISTA(maxloc(F_specBIS2(:,K,it),dim=1, mask=(.NOT.ISNAN(F_specBIS2(:,K,it)))))
      END IF
     ! WRITE(6,*) F_spec_tot(K,it)
    !  WRITE(6,*) Type(F_spec_tot)
    END DO ! DO K=1,NK

    DO AZ=1,NSTAAZI
        Amp_t_az_sum1(AZ,it)=SUM(F_spec1(AZ,:,it)*df,dim=1, mask =(.NOT.ISNAN (F_spec1(AZ,:,it)*df)))!just 1 dimension
        Amp_t_az_sum2(AZ,it)=SUM(F_spec2(AZ,:,it)*df,dim=1, mask = (.NOT.ISNAN (F_spec2(AZ,:,it)*df)))
        Amp_t_az_mf1(AZ,it)=maxval(F_spec1(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_spec1(AZ,:,it))))
        Amp_t_az_mf2(AZ,it)=maxval(F_spec2(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_spec2(AZ,:,it))))
        f_t_az_mamp1(AZ,it)=AZISTA(maxloc(F_spec1(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_spec1(AZ,:,it)))))
        f_t_az_mamp2(AZ,it)=AZISTA(maxloc(F_spec2(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_spec2(AZ,:,it)))))
        IF(FLAG_POWBIS)THEN
          Amp_t_az_sumBIS1(AZ,it)=SUM(F_specBIS1(AZ,:,it)*df,dim=1, mask = (.NOT.ISNAN (F_specBIS1(AZ,:,it)*df)))
          Amp_t_az_sumBIS2(AZ,it)=SUM(F_specBIS2(AZ,:,it)*df,dim=1, mask = (.NOT.ISNAN (F_specBIS2(AZ,:,it)*df)))
          Amp_t_az_mfBIS1(AZ,it)=maxval(F_specBIS1(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_specBIS1(AZ,:,it))))
          Amp_t_az_mfBIS2(AZ,it)=maxval(F_specBIS2(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_specBIS2(AZ,:,it))))
          f_t_az_mampBIS1(AZ,it)=AZISTA(maxloc(F_specBIS1(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_specBIS1(AZ,:,it)))))
          f_t_az_mampBIS2(AZ,it)=AZISTA(maxloc(F_specBIS2(AZ,:,it),dim=1, mask = (.NOT.ISNAN (F_specBIS2(AZ,:,it)))))
        END IF
      END DO
    DO Nkk=1,5 
      limf1 = 0.1*Nkk
      limf2 = 0.1*(Nkk+1)
      indf1=maxloc(freq1,mask=(freq1.LE.limf1),dim=1)!.AND.(freq2.GT.limf1)))
      indf2=maxloc(freq2,mask=(freq1.LE.limf2),dim=1)!.AND.(freq2.GT.limf2)))
      IF (indf1*indf2.NE.0) THEN
        IF(ALLOCATED(dfind)) DEALLOCATE(dfind)
        ALLOCATE(DFIND(indf2-indf1+1))
        IF (indf1.NE.indf2) THEN
          dfind=df(indf1:indf2)
          dfind(1)=freq2(indf1)-limf1
          dfind(indf2-indf1+1)=limf2-freq1(indf2)
        ELSE
          dfind=limf2-limf1
        END IF
      ELSE
        WRITE(6,*) 'Pb with indf1 or indf2'
      END IF
      DO AZ=1,NSTAAZI
        IF(indf1*indf2.NE.0) THEN
          IF (indf1.NE.indf2) THEN
            Amp_t_az_Delf1(AZ,Nkk,it)=SUM(F_spec1(AZ,indf1:indf2,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_spec1(AZ,indf1:indf2,it)*dfind)))
            Amp_t_az_Delf2(AZ,Nkk,it)=SUM(F_spec2(AZ,indf1:indf2,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_spec2(AZ,indf1:indf2,it)*dfind)))
            IF(FLAG_POWBIS)THEN
              Amp_t_az_DelfBIS1(AZ,Nkk,it)=SUM(F_specBIS1(AZ,indf1:indf2,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_specBIS1(AZ,indf1:indf2,it)*dfind)))
              Amp_t_az_DelfBIS2(AZ,Nkk,it)=SUM(F_specBIS2(AZ,indf1:indf2,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_specBIS2(AZ,indf1:indf2,it)*dfind)))
            END IF
          ELSE
            Amp_t_az_Delf1(AZ,Nkk,it)=SUM(F_spec1(AZ,indf1,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_spec1(AZ,indf1,it)*dfind)))
            Amp_t_az_Delf2(AZ,Nkk,it)=SUM(F_spec2(AZ,indf1,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_spec2(AZ,indf1,it)*dfind)))
            IF(FLAG_POWBIS)THEN
              Amp_t_az_DelfBIS1(AZ,Nkk,it)=SUM(F_specBIS1(AZ,indf1,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_specBIS1(AZ,indf1,it)*dfind)))
              Amp_t_az_DelfBIS2(AZ,Nkk,it)=SUM(F_specBIS2(AZ,indf1,it)*dfind,dim=1, mask = (.NOT.ISNAN (F_specBIS2(AZ,indf1,it)*dfind)))
            END IF
          END IF
        ELSE
          WRITE(6,*) 'Pb with indf1 or indf2'
        END IF
      END DO
    END DO

!-------------------------------------------------------------------------------
!---- Write for each time step -------------------------------------------------
! 
    iret = nf90_put_var(ncio1,varido1(11),Amp_t_f_sum1(:,it),start=(/1,IT/),count=(/NK,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_f_sum1')
    iret = nf90_put_var(ncio1,varido1(12),Amp_t_az_Delf1(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,5,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_az_Delf1')
    iret = nf90_put_var(ncio1,varido1(13),Amp_t_az_sum1(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_az_sum1')
    iret = nf90_put_var(ncio1,varido1(14),Amp_t_f_maz1(:,it),start=(/1,IT/), count=(/NK,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_f_maz1')
    iret = nf90_put_var(ncio1,varido1(15),Az_t_f_mamp1(:,it),start=(/1,IT/), count=(/NK,1/))
    call check_err(iret,'Error in nf90_put_var Az_t_f_mamp1')
    iret = nf90_put_var(ncio1,varido1(16),Amp_t_az_mf1(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_az_mf1')
    iret = nf90_put_var(ncio1,varido1(17),f_t_az_mamp1(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
    call check_err(iret,'Error in nf90_put_var f_t_az_mamp1')

    IF(FLAG_POWBIS)THEN
!      iret = nf90_put_var(ncio1,varido1(15),F_spec_totBIS1(:,it),start=(/1,IT/),count=(/NK,1/))
!      call check_err(iret,'Error in nf90_put_var F_spec_tot')
!      iret = nf90_put_var(ncio1,varido1(16),F_specBIS1(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,NK,1/))
!      call check_err(iret,'Error in nf90_put_var F_spec')
!      iret = nf90_put_var(ncio1,varido1(17),DISTMBIS1(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,NK,1/))
!      call check_err(iret,'Error in nf90_put_var DISTM')
!      iret = nf90_put_var(ncio1,varido1(18),AMPMBIS1(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,NK,1/))
!      call check_err(iret,'Error in nf90_put_var AMPM')
      iret = nf90_put_var(ncio1,varido1(18),Amp_t_f_sumBIS1(:,it),start=(/1,IT/),count=(/NK,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_f_sumBIS1')
      iret = nf90_put_var(ncio1,varido1(19),Amp_t_az_DelfBIS1(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,5,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_az_DelfBIS1')
      iret = nf90_put_var(ncio1,varido1(20),Amp_t_az_sumBIS1(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_az_sumBIS1')
      iret = nf90_put_var(ncio1,varido1(21),Amp_t_f_mazBIS1(:,it),start=(/1,IT/), count=(/NK,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_f_mazBIS1')
      iret = nf90_put_var(ncio1,varido1(22),Az_t_f_mampBIS1(:,it),start=(/1,IT/), count=(/NK,1/))
      call check_err(iret,'Error in nf90_put_var Az_t_f_mampBIS1')
      iret = nf90_put_var(ncio1,varido1(23),Amp_t_az_mfBIS1(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_az_mfBIS1')
      iret = nf90_put_var(ncio1,varido1(24),f_t_az_mampBIS1(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
      call check_err(iret,'Error in nf90_put_var f_t_az_mampBIS1')
    END IF
    ! Save for the PropaNoWind file -------------------------------------------------------------
    iret = nf90_put_var(ncio2,varido2(11),Amp_t_f_sum2(:,it),start=(/1,IT/),count=(/NK,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_f_sum2')
    iret = nf90_put_var(ncio2,varido2(12),Amp_t_az_Delf2(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,5,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_az_Delf2')
    iret = nf90_put_var(ncio2,varido2(13),Amp_t_az_sum2(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_az_sum2')
    iret = nf90_put_var(ncio2,varido2(14),Amp_t_f_maz2(:,it),start=(/1,IT/), count=(/NK,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_f_maz2')
    iret = nf90_put_var(ncio1,varido2(15),Az_t_f_mamp2(:,it),start=(/1,IT/), count=(/NK,1/))
    call check_err(iret,'Error in nf90_put_var Az_t_f_mamp2')
    iret = nf90_put_var(ncio1,varido2(16),Amp_t_az_mf2(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
    call check_err(iret,'Error in nf90_put_var Amp_t_az_mf2')
    iret = nf90_put_var(ncio1,varido2(17),f_t_az_mamp2(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
    call check_err(iret,'Error in nf90_put_var f_t_az_mamp2')

    IF(FLAG_POWBIS)THEN
!      iret = nf90_put_var(ncio2,varido2(15),F_spec_totBIS2(:,it),start=(/1,IT/),count=(/NK,1/))
!      call check_err(iret,'Error in nf90_put_var F_spec_tot')
!      iret = nf90_put_var(ncio2,varido2(16),F_specBIS2(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,NK,1/))
!      call check_err(iret,'Error in nf90_put_var F_spec')
!      iret = nf90_put_var(ncio2,varido2(17),DISTMBIS2(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,NK,1/))
!      call check_err(iret,'Error in nf90_put_var DISTM')
!      iret = nf90_put_var(ncio2,varido2(18),AMPMBIS2(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,NK,1/))
!      call check_err(iret,'Error in nf90_put_var AMPM')
      iret = nf90_put_var(ncio2,varido2(18),Amp_t_f_sumBIS2(:,it),start=(/1,IT/),count=(/NK,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_f_sumBIS2')
      iret = nf90_put_var(ncio2,varido2(19),Amp_t_az_DelfBIS2(:,:,it),start=(/1,1,IT/), count=(/NSTAAZI,5,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_az_DelfBIS2')
      iret = nf90_put_var(ncio2,varido2(20),Amp_t_az_sumBIS2(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_az_sumBIS2')
      iret = nf90_put_var(ncio2,varido2(21),Amp_t_f_mazBIS2(:,it),start=(/1,IT/), count=(/NK,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_f_mazBIS2')
      iret = nf90_put_var(ncio2,varido2(22),Az_t_f_mampBIS2(:,it),start=(/1,IT/), count=(/NK,1/))
      call check_err(iret,'Error in nf90_put_var Az_t_f_mampBIS2')
      iret = nf90_put_var(ncio2,varido2(23),Amp_t_az_mfBIS2(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
      call check_err(iret,'Error in nf90_put_var Amp_t_az_mfBIS2')
      iret = nf90_put_var(ncio2,varido2(24),f_t_az_mampBIS2(:,it),start=(/1,IT/), count=(/NSTAAZI,1/))
      call check_err(iret,'Error in nf90_put_var f_t_az_mampBIS2')
    END IF

  END DO ! DO it=1,NT
  iret = nf90_put_var(ncio1,varido1(1),times)
  call check_err(iret,'Error in nf90_put_var times')
  iret = nf90_close(ncio1)
  call check_err(iret,'Error in closing output file')
  iret = nf90_put_var(ncio2,varido2(1),times)
  call check_err(iret,'Error in nf90_put_var times')
  iret = nf90_close(ncio2)
  call check_err(iret,'Error in closing output file')
END DO ! DO I=1, NSTATIONS

WRITE(NDSO,*) ' *** END OF PROGRAM *** '
END PROGRAM

!===============================================================================
!-------------------------------------------------------------------------------
  SUBROUTINE attenuation(DIS,f,Vratio,ATTTMP,LL)
!      Given a distance from the source 'DIS' (in km), a frequency 
!  'f' (in Hz) and a wind value at the station (considering wind is 
!  constant all along the path) 'Vratio', this function will 
!  calculate the Attenuation coefficient of pressure wave 'ATTTMP'
   use netcdf
   IMPLICIT NONE

   INTEGER                         :: LL, f1, F2
   REAL, DIMENSION(LL)             :: Vratio, beta_val
   DOUBLE PRECISION, DIMENSION(LL) :: DIS
   REAL                            :: f_table(6), Veff_table(12), alpha_table(6),sigma_table(6), beta_val0(12)
   REAL                            :: beta_table(12,6)
   REAL, DIMENSION(LL,12)          :: Vratio2,Veff_table2
   INTEGER, DIMENSION(LL)          :: V1,V2
   DOUBLE PRECISION, ALLOCATABLE   :: ATTTMP(:)
   REAL                            :: f, delta, MinVeff, MaxVeff, alpha_val,sigma_val

    IF (ALLOCATED(ATTTMP)) DEALLOCATE(ATTTMP)
    ALLOCATE(ATTTMP(LL))
    !WRITE(6,*) Vratio
    f_table = (/ 0.1, 0.2, 0.4, 0.8, 1.6, 3.2 /)
    Veff_table = (/ 0.85, 0.88, 0.91, 0.94, 0.97, 1.0, 1.03, 1.06, 1.09, 1.12, 1.15, 1.18 /)
    MinVeff = 0.85
    MaxVeff = 1.18
    WHERE (Vratio.LT.MinVeff) Vratio=MinVeff
    WHERE (Vratio.GT.MaxVeff) Vratio=MaxVeff

! Table for alpha
    alpha_table = (/ -0.28, -0.33, -0.39, -0.47, -0.59, -0.69 /)
! Table for beta
    beta_table(1,:) = (/ -1.0, -1.2, -1.35, -1.70, -1.95, -2.4 /)
    beta_table(2,:) = (/ -1.0, -1.2, -1.4, -1.6, -1.95, -2.30 /)
    beta_table(3,:) = (/ -1.0, -1.15, -1.4, -1.6, -1.85, -2.1 /)
    beta_table(4,:) = (/ -1.0, -1.1, -1.35, -1.55, -1.75, -1.85 /)
    beta_table(5,:) = (/ -1.0, -1.05, -1.15, -1.30, -1.40, -1.45 /)
    beta_table(6,:) = (/ -0.90, -0.95, -1.0, -1.05, -1.15, -1.20 /)
    beta_table(7,:) = (/ -0.90, -0.90, -0.90, -0.90, -0.95, -1.0 /)
    beta_table(8,:) = (/ -0.85, -0.85, -0.90, -0.90, -0.95, -1.0 /)
    beta_table(9,:) = (/ -0.85, -0.85, -0.90, -0.90, -1.0, -1.05 /)
    beta_table(10,:) = (/ -0.85, -0.85, -0.85, -0.90, -1.0, -1.10 /)
    beta_table(11,:) = (/ -0.80, -0.80, -0.85, -0.95, -1.05, -1.15 /)
    beta_table(12,:) = (/ -0.80, -0.80, -0.85, -0.90, -1.05, -1.20 /)
! Table for sigma
    sigma_table = (/ 79, 55, 43, 36, 27, 20 /);
! Value for delta
    delta = 180 ;

! Interp for f values
! .....................
    f1=maxloc(f_table,dim=1,mask=(f_table<f))
    f2=minloc(f_table,dim=1,mask=(f_table>f))
    If(f1.EQ.0) THEN
! This is f<f_table and you do 'expand' with same function as for the border
      f1=1 !f_table(1)
      f2=2 !f_table(2)
    ELSEIF (f2.EQ.0)THEN
! This is f>f_table and you do 'expand' with same function as for the border
      f1=5 ! f_table(5)
      f2=6 ! f_table(6)
    END IF

    alpha_val = ((f_table(f2)*alpha_table(f1) - f_table(f1)*alpha_table(f2))/(f_table(f2) &
     -f_table(f1))) + (f*(alpha_table(f2) - alpha_table(f1))/(f_table(f2)-f_table(f1)))

    sigma_val = (f_table(f2)*sigma_table(f1) - f_table(f1)*sigma_table(f2))/(f_table(f2) &
     -f_table(f1)) + f*(sigma_table(f2) - sigma_table(f1))/(f_table(f2)-f_table(f1))

    beta_val0(:) = (f_table(f2)*beta_table(:,f1) - f_table(f1)*beta_table(:,f2))/ & 
(f_table(f2)-f_table(f1)) + f*(beta_table(:,f2) - beta_table(:,f1))/(f_table(f2)-f_table(f1))

! Interp for Veff values
! .................... 
! WARNING CHECK TEST
    Veff_table2 = SPREAD(source=Veff_table, dim=1,ncopies=LL)
    Vratio2 = SPREAD(source=Vratio,dim=2,ncopies=12)

    V1=maxloc(Veff_table2,mask=Veff_table2<Vratio2,dim=2)
    V2=minloc(Veff_table2,mask=Veff_table2>Vratio2,dim=2)

!    IF(V1.EQ.0) V2=2
!    IF(V1.EQ.0) V1=1
    WHERE (V1.EQ.0) V2=2
    WHERE (V1.EQ.0) V1=1
! This is V<V_table and you do 'expand' with same function as for the border
!    IF(V2.EQ.0) V1=11
!    IF(V2.EQ.0) V2=12
    WHERE (V2.EQ.0) V1=11
    WHERE (V2.EQ.0) V2=12
! This is f>f_table and you do 'expand' with same function as for the border

    beta_val(:) = ((Veff_table(V2)*beta_val0(V1) - Veff_table(V1)*beta_val0(V2))/(Veff_table(V2) &
      -Veff_table(V1)) + Vratio(:)*(beta_val0(V2) - beta_val0(V1))/(Veff_table(V2)-Veff_table(V1)))

    ATTTMP(:) = ((1.)/DIS(:))*(10)**((alpha_val)*DIS(:)/(20)) + &
     (DIS(:)**(beta_val(:)))/((1)+(10)**(((delta) - DIS(:))/(sigma_val)))
   return

  END SUBROUTINE
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
  SUBROUTINE attenuationAZE (DIS,f,Vraze,ATTTMP,LL)
!      Given a distance from the source 'DIS' (in km), a frequency 
!  'f' (in Hz) and a wind value at the station (considering wind is 
!  constant all along the path) 'Vratioaz', this function will 
!  calculate the Attenuation coefficient of pressure wave 'ATTTMP'


   IMPLICIT NONE

   INTEGER                         :: LL, f1, F2, V1,V2
   DOUBLE PRECISION, DIMENSION(LL) :: DIS
   REAL                            :: f_table(6), Veff_table(12), alpha_table(6),sigma_table(6), beta_val0(12),Vraze
   REAL                            :: beta_table(12,6)
   DOUBLE PRECISION, ALLOCATABLE   :: ATTTMP(:)
   REAL                            :: f, delta, MinVeff, MaxVeff, alpha_val, beta_val, sigma_val

    IF (ALLOCATED(ATTTMP)) DEALLOCATE(ATTTMP)
    ALLOCATE(ATTTMP(LL))
  
    f_table = (/ 0.1, 0.2, 0.4, 0.8, 1.6, 3.2 /)
    Veff_table = (/ 0.85, 0.88, 0.91, 0.94, 0.97, 1.0, 1.03, 1.06, 1.09, 1.12, 1.15, 1.18 /)
    MinVeff = minval(Veff_table)
    MaxVeff = maxval(Veff_table)

    IF(Vraze.LT.MinVeff) Vraze=MinVeff
    IF(Vraze.GT.MaxVeff) Vraze=MaxVeff

! Table for alpha
    alpha_table = (/ -0.28, -0.33, -0.39, -0.47, -0.59, -0.69 /)
! Table for beta
    beta_table(1,:) = (/ -1.0, -1.2, -1.35, -1.70, -1.95, -2.4 /)
    beta_table(2,:) = (/ -1.0, -1.2, -1.4, -1.6, -1.95, -2.30 /)
    beta_table(3,:) = (/ -1.0, -1.15, -1.4, -1.6, -1.85, -2.1 /)
    beta_table(4,:) = (/ -1.0, -1.1, -1.35, -1.55, -1.75, -1.85 /)
    beta_table(5,:) = (/ -1.0, -1.05, -1.15, -1.30, -1.40, -1.45 /)
    beta_table(6,:) = (/ -0.90, -0.95, -1.0, -1.05, -1.15, -1.20 /)
    beta_table(7,:) = (/ -0.90, -0.90, -0.90, -0.90, -0.95, -1.0 /)
    beta_table(8,:) = (/ -0.85, -0.85, -0.90, -0.90, -0.95, -1.0 /)
    beta_table(9,:) = (/ -0.85, -0.85, -0.90, -0.90, -1.0, -1.05 /)
    beta_table(10,:) = (/ -0.85, -0.85, -0.85, -0.90, -1.0, -1.10 /)
    beta_table(11,:) = (/ -0.80, -0.80, -0.85, -0.95, -1.05, -1.15 /)
    beta_table(12,:) = (/ -0.80, -0.80, -0.85, -0.90, -1.05, -1.20 /)
! Table for sigma
    sigma_table = (/ 79, 55, 43, 36, 27, 20 /);
! Value for delta
    delta = 180 ;

! Interp for f values
! .....................
    f1=maxloc(f_table,dim=1,mask=(f_table<f))
    f2=minloc(f_table,dim=1,mask=(f_table>f))
    If(f1.EQ.0) THEN
! This is f<f_table and you do 'expand' with same function as for the border
      f1=1 !f_table(1)
      f2=2 !f_table(2)
    ELSEIF (f2.EQ.0)THEN
! This is f>f_table and you do 'expand' with same function as for the border
      f1=5 ! f_table(5)
      f2=6 ! f_table(6)
    END IF

    alpha_val = (f_table(f2)*alpha_table(f1) - f_table(f1)*alpha_table(f2))/(f_table(f2)-f_table(f1)) + &
 f*(alpha_table(f2) - alpha_table(f1))/(f_table(f2)-f_table(f1))

    sigma_val = (f_table(f2)*sigma_table(f1) - f_table(f1)*sigma_table(f2))/(f_table(f2)-f_table(f1)) + &
 f*(sigma_table(f2) - sigma_table(f1))/(f_table(f2)-f_table(f1))

    beta_val0(:) = (f_table(f2)*beta_table(:,f1) - f_table(f1)*beta_table(:,f2))/(f_table(f2)-f_table(f1)) +&
 f*(beta_table(:,f2) - beta_table(:,f1))/(f_table(f2)-f_table(f1))

! Interp for Veff values
! .................... 

    V1=maxloc(Veff_table,dim=1,mask=Veff_table<Vraze)
    V2=minloc(Veff_table,dim=1,mask=Veff_table>Vraze)
    IF(V1.EQ.0) V2=2
    IF(V1.EQ.0) V1=1
! This is f<f_table and you do 'expand' with same function as for the border
    IF(V2.EQ.0) V1=11
    IF(V2.EQ.0) V2=12
! This is f>f_table and you do 'expand' with same function as for the border

    beta_val = ((Veff_table(V2)*beta_val0(V1) - Veff_table(V1)*beta_val0(V2))/(Veff_table(V2) &
     -Veff_table(V1)) + Vraze*(beta_val0(V2) - beta_val0(V1))/(Veff_table(V2)-Veff_table(V1)))

    ATTTMP(:) = ((1.)/DIS(:))*(10)**((alpha_val)*DIS(:)/(20)) + (DIS(:)**(beta_val))&
/((1)+(10)**(((delta) - DIS(:))/(sigma_val)));

   return

  END SUBROUTINE
!-------------------------------------------------------------------------------
!===============================================================================
!-------------------------------------------------------------------------------
  SUBROUTINE attenuationNoWind(DIS,f,ATTTMP,LL)
!      Given a distance from the source 'DIS' (in km), a frequency 
!  'f' (in Hz) this function will 
!  calculate the Attenuation coefficient of pressure wave 'ATTTMP'(the wind ratio is considered to be 1)
   use netcdf
   IMPLICIT NONE

   INTEGER                         :: LL, f1, F2
   DOUBLE PRECISION, DIMENSION(LL) :: DIS
   REAL                            :: f_table(6), alpha_table(6),sigma_table(6)
   REAL                            :: beta_table(6)
   DOUBLE PRECISION, ALLOCATABLE   :: ATTTMP(:)
   REAL                            :: f, delta, alpha_val,sigma_val, beta_val

    IF (ALLOCATED(ATTTMP)) DEALLOCATE(ATTTMP)
    ALLOCATE(ATTTMP(LL))
    !WRITE(6,*) Vratio
    f_table = (/ 0.1, 0.2, 0.4, 0.8, 1.6, 3.2 /)
! Table for alpha
    alpha_table = (/ -0.28, -0.33, -0.39, -0.47, -0.59, -0.69 /)
! Table for beta
    beta_table = (/ -0.90, -0.95, -1.0, -1.05, -1.15, -1.20 /)
! Table for sigma
    sigma_table = (/ 79, 55, 43, 36, 27, 20 /);
! Value for delta
    delta = 180 ;

! Interp for f values
! .....................
    f1=maxloc(f_table,dim=1,mask=(f_table<f))
    f2=minloc(f_table,dim=1,mask=(f_table>f))
    If(f1.EQ.0) THEN
! This is f<f_table and you do 'expand' with same function as for the border
      f1=1 !f_table(1)
      f2=2 !f_table(2)
    ELSEIF (f2.EQ.0)THEN
! This is f>f_table and you do 'expand' with same function as for the border
      f1=5 ! f_table(5)
      f2=6 ! f_table(6)
    END IF

    alpha_val = ((f_table(f2)*alpha_table(f1) - f_table(f1)*alpha_table(f2))/(f_table(f2) &
     -f_table(f1))) + (f*(alpha_table(f2) - alpha_table(f1))/(f_table(f2)-f_table(f1)))

    sigma_val = (f_table(f2)*sigma_table(f1) - f_table(f1)*sigma_table(f2))/(f_table(f2) &
     -f_table(f1)) + f*(sigma_table(f2) - sigma_table(f1))/(f_table(f2)-f_table(f1))

    beta_val = (f_table(f2)*beta_table(f1) - f_table(f1)*beta_table(f2))/ & 
(f_table(f2)-f_table(f1)) + f*(beta_table(f2) - beta_table(f1))/(f_table(f2)-f_table(f1))

    ATTTMP(:) = ((1.)/DIS(:))*(10)**((alpha_val)*DIS(:)/(20)) + &
     (DIS(:)**(beta_val))/((1)+(10)**(((delta) - DIS(:))/(sigma_val)))
   return

  END SUBROUTINE
!-------------------------------------------------------------------------------
!===============================================================================!-------------------------------------------------------------------------------
SUBROUTINE DISTAZ (la1,lo1,la2,lo2,az2from1,az1from2)
! Given 2 points : 1(la1,lo1) and 2(la2,lo2) we compute two angles: az2from1 : azimuth of point 2 seen from point 1, az1from2: azimuth of point 1 seen from point2
implicit none
REAL, INTENT(IN)          :: la1(1),lo1(1),la2,lo2
REAL, INTENT(OUT)         :: az2from1, az1from2
REAL                      :: sa,ssa,sca,so,ea,eo,sea,eca,ecs,ses,rad

rad=4*datan(1d0)/180.

sa=atan(.99327*tan(la1(1)*rad));
ssa=sin(sa);
sca=cos(sa);
so=lo1(1)*rad;
ea=atan(.99327*tan(la2*rad));
eo=lo2*rad;
sea=sin(ea);
eca=cos(ea);
ecs=cos(eo-so);
ses=sin(eo-so);
      
az1from2=atan2(-ses,((ssa/sca*eca)-(sea*ecs)));
az2from1=atan2(ses,((sea*sca/eca)-(ssa*ecs)));
      
az1from2=az1from2/rad;
az2from1=az2from1/rad;
      
IF (az1from2.LT.0) az1from2 =az1from2 +360

IF (az2from1.LT.0) az2from1 =az2from1 +360

return
END SUBROUTINE
!-------------------------------------------------------------------------------
!===============================================================================
!---------------------------------------------------------------------
SUBROUTINE W3CRNC (NCTYPE,ncfile, ncid,dimid1, dimln1, varid, FLAG_bis)
  USE W3GDATMD
  USE NETCDF

  implicit none
  
  CHARACTER*(256), INTENT(IN)    :: ncfile
  INTEGER, INTENT(IN)            :: dimln1(6), NCTYPE
  INTEGER, INTENT(OUT)           :: dimid1(6), varid(28), ncid
  INTEGER                        :: iret
  LOGICAL                        :: FLAG_bis

! 
! Creation in NetCDF3 or netCDF4
!  
  IF (NCTYPE.EQ.3) iret = nf90_create(trim(ncfile), NF90_CLOBBER, ncid)
  IF (NCTYPE.EQ.4) iret = nf90_create(trim(ncfile), NF90_NETCDF4, ncid)

! 
! Define dimensions
! 
  iret = nf90_def_dim(ncid, 'time', dimln1(1), dimid1(1))
  iret = nf90_def_dim(ncid, 'station', dimln1(2), dimid1(2))
  iret = nf90_def_dim(ncid, 'string16', dimln1(3), dimid1(3))
  iret = nf90_def_dim(ncid, 'frequency', dimln1(4), dimid1(4))
  iret = nf90_def_dim(ncid, 'direction', dimln1(5), dimid1(5))
  iret = nf90_def_dim(ncid, 'freq_band', dimln1(6), dimid1(6))
! 
! Define variables
! 

! Time 
  iret = nf90_def_var(ncid, 'time', NF90_DOUBLE, dimid1(1), varid(1))
  
  iret = nf90_put_att(ncid, varid(1),'long_name','julian day (UT)')
  iret = nf90_put_att(ncid, varid(1),'standard_name','time')
  iret = nf90_put_att(ncid, varid(1),'units','days since 1990-01-01 00:00:00')
  iret = nf90_put_att(ncid, varid(1),'conventions','Relative julian days with decimal part (as parts of the day)')
  iret = nf90_put_att(ncid, varid(1),'axis','T')

! station_name
  iret = nf90_def_var(ncid, 'station_name', NF90_CHAR, (/dimid1(3),dimid1(2)/), varid(2))

  iret = nf90_put_att(ncid, varid(2),'long_name','station name')
  iret = nf90_put_att(ncid, varid(2),'content', 'XW')
  iret = nf90_put_att(ncid, varid(2),'associates','station string16')

! station
  iret = nf90_def_var(ncid, 'station', NF90_INT, (/dimid1(2)/), varid(27))

  iret = nf90_put_att(ncid, varid(27),'long_name','station id')
  iret = nf90_put_att(ncid, varid(27),'_FillValue',NF90_FILL_INT)
  iret = nf90_put_att(ncid, varid(27),'axis','X')

! string16
  iret = nf90_def_var(ncid, 'string16', NF90_INT, (/dimid1(3)/), varid(28))

  iret = nf90_put_att(ncid, varid(28),'long_name','station_name number of characters')
  iret = nf90_put_att(ncid, varid(28),'_FillValue',NF90_FILL_INT)
  iret = nf90_put_att(ncid, varid(28),'axis','W')

! longitude
  iret = nf90_def_var(ncid, 'longitude', NF90_FLOAT, (/dimid1(2),dimid1(1)/),varid(4))

  iret = nf90_put_att(ncid,varid(4),'long_name','longitude')
  iret = nf90_put_att(ncid,varid(4),'standard_name','longitude')
  iret = nf90_put_att(ncid,varid(4),'globwave_name','longitude')
  iret = nf90_put_att(ncid,varid(4),'units','degree_east')
  iret = nf90_put_att(ncid,varid(4),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(4),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(4),'valid_min',-180.0)
  iret = nf90_put_att(ncid,varid(4),'valid_max',360.0)
  iret = nf90_put_att(ncid,varid(4),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(4),'content','TX')
  iret = nf90_put_att(ncid,varid(4),'associates','time station')

! latitude
  iret = nf90_def_var(ncid, 'latitude', NF90_FLOAT, (/dimid1(2),dimid1(1)/),varid(5))

  iret = nf90_put_att(ncid,varid(5),'long_name','latitude')
  iret = nf90_put_att(ncid,varid(5),'standard_name','latitude')
  iret = nf90_put_att(ncid,varid(5),'globwave_name','latitude')
  iret = nf90_put_att(ncid,varid(5),'units','degree_north')
  iret = nf90_put_att(ncid,varid(5),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(5),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(5),'valid_min',-90.0)
  iret = nf90_put_att(ncid,varid(5),'valid_max',180.0)
  iret = nf90_put_att(ncid,varid(5),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(5),'content','TX')
  iret = nf90_put_att(ncid,varid(5),'associates','time station')

! 
! Define variables
! 

! frequency
  iret = nf90_def_var(ncid, 'frequency', NF90_FLOAT,(/dimid1(4)/), varid(6))
  call check_err(iret, 'Error in nf90_def_var frequency')
  iret = nf90_put_att(ncid,varid(6),'long_name','frequency of center band')
  iret = nf90_put_att(ncid,varid(6),'standard_name','sea_surface_wave_frequency')
  iret = nf90_put_att(ncid,varid(6),'globwave_name','frequency')
  iret = nf90_put_att(ncid,varid(6),'units','s-1')
  iret = nf90_put_att(ncid,varid(6),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(6),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(6),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(6),'valid_max',10.)
  iret = nf90_put_att(ncid,varid(6),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(6),'axis','Y')

! frequency band
  iret = nf90_def_var(ncid, 'freq_band', NF90_FLOAT,(/dimid1(6)/), varid(26))
  call check_err(iret, 'Error in nf90_def_var frequency band')
  iret = nf90_put_att(ncid,varid(26),'long_name','frequency band for summation (lower frequency, delt = 0.1Hz)')
  iret = nf90_put_att(ncid,varid(26),'standard_name','frequency band for summation')
  iret = nf90_put_att(ncid,varid(26),'globwave_name','frequency band')
  iret = nf90_put_att(ncid,varid(26),'units','s-1')
  iret = nf90_put_att(ncid,varid(26),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(26),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(26),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(26),'valid_max',10.)
  iret = nf90_put_att(ncid,varid(26),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(26),'axis','Z')
! frequency 1
  iret = nf90_def_var(ncid, 'frequency1', NF90_FLOAT,(/dimid1(4)/), varid(7))

  iret = nf90_put_att(ncid,varid(7),'long_name','frequency of lower band')
  iret = nf90_put_att(ncid,varid(7),'standard_name','frequency_of_lower_band')
  iret = nf90_put_att(ncid,varid(7),'globwave_name','frequency_lower_band')
  iret = nf90_put_att(ncid,varid(7),'units','s-1')
  iret = nf90_put_att(ncid,varid(7),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(7),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(7),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(7),'valid_max',10.)
  iret = nf90_put_att(ncid,varid(7),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(7),'content','Y')
  iret = nf90_put_att(ncid,varid(7),'associates','frequency')
! frequency 2
  iret = nf90_def_var(ncid, 'frequency2', NF90_FLOAT,(/dimid1(4)/), varid(8))

  iret = nf90_put_att(ncid,varid(8),'long_name','frequency of upper band')
  iret = nf90_put_att(ncid,varid(8),'standard_name','frequency_of_upper_band')
  iret = nf90_put_att(ncid,varid(8),'globwave_name','frequency_upper_band')
  iret = nf90_put_att(ncid,varid(8),'units','s-1')
  iret = nf90_put_att(ncid,varid(8),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(8),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(8),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(8),'valid_max',10.)
  iret = nf90_put_att(ncid,varid(8),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(8),'content','Y')
  iret = nf90_put_att(ncid,varid(8),'associates','frequency')

! direction
  iret = nf90_def_var(ncid, 'direction',NF90_FLOAT,(/dimid1(5)/),varid(9))
  call check_err(iret,'Errorin nf90_def_var direction')
  iret = nf90_put_att(ncid,varid(9),'long_name','back_azimuth')
  iret = nf90_put_att(ncid,varid(9),'standard_name','azimuth')
  iret = nf90_put_att(ncid,varid(9),'globwave_name','azimuth')
  iret = nf90_put_att(ncid,varid(9),'units','deg')
  iret = nf90_put_att(ncid,varid(9),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(9),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(9),'valid_min',-180.0)
  iret = nf90_put_att(ncid,varid(9),'valid_max',360.0)
  iret = nf90_put_att(ncid,varid(9),'_FillValue',NF90_FILL_FLOAT)
  iret = nf90_put_att(ncid,varid(9),'axis','A')

!edzf
  iret = nf90_def_var(ncid, 'edzf',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(11))

  iret = nf90_put_att(ncid,varid(11),'long_name','ground acoustic power freq spectral density - sum over azi of F_spec')
  iret = nf90_put_att(ncid,varid(11),'standard_name','ground acoustic power freq spectral density')
  iret = nf90_put_att(ncid,varid(11),'globwave_name','ground acoustic power freq spectral density')
  iret = nf90_put_att(ncid,varid(11),'units','m2 s')
  iret = nf90_put_att(ncid,varid(11),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(11),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(11),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(11),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(11),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(11),'content','TY')
  iret = nf90_put_att(ncid,varid(11),'associates','time frequency')

!edzfth
  iret = nf90_def_var(ncid, 'edzfth',NF90_DOUBLE,(/dimid1(5),dimid1(6),dimid1(1)/),varid(12))

  iret = nf90_put_att(ncid,varid(12),'long_name','ground acoustic directional power spectral density - sum over freq of F_spec for frequency bands ')
  iret = nf90_put_att(ncid,varid(12),'standard_name','ground acoustic directional power spectral density (freq_band)')
  iret = nf90_put_att(ncid,varid(12),'globwave_name','ground acoustic directional power spectral density (freq_band)')
  iret = nf90_put_att(ncid,varid(12),'units','m2 s')
  iret = nf90_put_att(ncid,varid(12),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(12),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(12),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(12),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(12),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(12),'content','TZA')
  iret = nf90_put_att(ncid,varid(12),'associates','direction time frequency')

!edzth
  iret = nf90_def_var(ncid, 'edzth',NF90_DOUBLE,(/dimid1(5),dimid1(1)/),varid(13))

  iret = nf90_put_att(ncid,varid(13),'long_name','ground acoustic power azi spectral density - sum over freq of F_spec')
  iret = nf90_put_att(ncid,varid(13),'standard_name','ground acoustic power azi spectral density')
  iret = nf90_put_att(ncid,varid(13),'globwave_name','ground acoustic power azi spectral density')
  iret = nf90_put_att(ncid,varid(13),'units','m2 s')
  iret = nf90_put_att(ncid,varid(13),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(13),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(13),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(13),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(13),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(13),'content','TA')
  iret = nf90_put_att(ncid,varid(13),'associates','direction time')

!'edzmaz' edz for max of F_spec over AZ
  iret = nf90_def_var(ncid, 'edzmaz',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(14))

  iret = nf90_put_att(ncid,varid(14),'long_name','max over az of F_spec (ground acoustic power spectral density)')
  iret = nf90_put_att(ncid,varid(14),'standard_name','max of F_spec over az')
  iret = nf90_put_att(ncid,varid(14),'globwave_name','max of F_spec over az')
  iret = nf90_put_att(ncid,varid(14),'units','m2 s')
  iret = nf90_put_att(ncid,varid(14),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(14),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(14),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(14),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(14),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(14),'content','TY')
  iret = nf90_put_att(ncid,varid(14),'associates','frequency time')
!'azmax' az corres. to max of F_spec over Az
  iret = nf90_def_var(ncid, 'azmax',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(15))

  iret = nf90_put_att(ncid,varid(15),'long_name','azi of max - max of F_spec over freq')
  iret = nf90_put_att(ncid,varid(15),'standard_name','azi of max')
  iret = nf90_put_att(ncid,varid(15),'globwave_name','azi of max')
  iret = nf90_put_att(ncid,varid(15),'units','m2 s')
  iret = nf90_put_att(ncid,varid(15),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(15),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(15),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(15),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(15),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(15),'content','TY')
  iret = nf90_put_att(ncid,varid(15),'associates','frequency time')
!'edzmf' edz for max of F_spec over f
  iret = nf90_def_var(ncid, 'edzmf',NF90_DOUBLE,(/dimid1(5),dimid1(1)/),varid(16))

  iret = nf90_put_att(ncid,varid(16),'long_name','max over freq of F_spec (ground acoustic power spectral density)')
  iret = nf90_put_att(ncid,varid(16),'standard_name','max of F_spec over freq')
  iret = nf90_put_att(ncid,varid(16),'globwave_name','max of F_spec over freq')
  iret = nf90_put_att(ncid,varid(16),'units','m2 s')
  iret = nf90_put_att(ncid,varid(16),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(16),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(16),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(16),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(16),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(16),'content','TA')
  iret = nf90_put_att(ncid,varid(16),'associates','direction time')
!'fmax' freq corres. to max of F_spec over f
  iret = nf90_def_var(ncid, 'fmax',NF90_DOUBLE,(/dimid1(5),dimid1(1)/),varid(17))

  iret = nf90_put_att(ncid,varid(17),'long_name','freq of max - max of F_spec over freq')
  iret = nf90_put_att(ncid,varid(17),'standard_name','freq of max')
  iret = nf90_put_att(ncid,varid(17),'globwave_name','freq of max')
  iret = nf90_put_att(ncid,varid(17),'units','m2 s')
  iret = nf90_put_att(ncid,varid(17),'scale_factor',1.)
  iret = nf90_put_att(ncid,varid(17),'add_offset',0.)
  iret = nf90_put_att(ncid,varid(17),'valid_min',0.)
  iret = nf90_put_att(ncid,varid(17),'valid_max',1.E20)
  iret = nf90_put_att(ncid,varid(17),'_FillValue',NF90_FILL_DOUBLE)
  iret = nf90_put_att(ncid,varid(17),'content','TA')
  iret = nf90_put_att(ncid,varid(17),'associates','direction time')

!OLD OUTPUTS
!dist_max
!  iret = nf90_def_var(ncid, 'dist_source_max',NF90_DOUBLE,(/dimid1(5),dimid1(4),dimid1(1)/),varid(13))

!  iret = nf90_put_att(ncid,varid(13),'long_name','distance of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(13),'standard_name','distance of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(13),'globwave_name','distance of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(13),'units','m2 s')
!  iret = nf90_put_att(ncid,varid(13),'scale_factor',1.)
!  iret = nf90_put_att(ncid,varid(13),'add_offset',0.)
!  iret = nf90_put_att(ncid,varid(13),'valid_min',0.)
!  iret = nf90_put_att(ncid,varid(13),'valid_max',1.E20)
!  iret = nf90_put_att(ncid,varid(13),'_FillValue',NF90_FILL_DOUBLE)
!  iret = nf90_put_att(ncid,varid(13),'content','TYA')
!  iret = nf90_put_att(ncid,varid(13),'associates','direction time frequency')

!amp_max
!  iret = nf90_def_var(ncid, 'source_max',NF90_DOUBLE,(/dimid1(5),dimid1(4),dimid1(1)/),varid(14))

!  iret = nf90_put_att(ncid,varid(14),'long_name','amplitude of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(14),'standard_name','amplitude of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(14),'globwave_name','amplitude of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(14),'units','m2 s')
!  iret = nf90_put_att(ncid,varid(14),'scale_factor',1.)
!  iret = nf90_put_att(ncid,varid(14),'add_offset',0.)
!  iret = nf90_put_att(ncid,varid(14),'valid_min',0.)
!  iret = nf90_put_att(ncid,varid(14),'valid_max',1.E20)
!  iret = nf90_put_att(ncid,varid(14),'_FillValue',NF90_FILL_DOUBLE)
!  iret = nf90_put_att(ncid,varid(14),'content','TYA')
!  iret = nf90_put_att(ncid,varid(14),'associates','direction time frequency')

  IF(FLAG_BIS)THEN
  !edzf
    iret = nf90_def_var(ncid, 'edzf_bis',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(18))

    iret = nf90_put_att(ncid,varid(18),'long_name','ground acoustic power freq spectral density - sum over azi of F_spec')
    iret = nf90_put_att(ncid,varid(18),'standard_name','ground acoustic power freq spectral density')
    iret = nf90_put_att(ncid,varid(18),'globwave_name','ground acoustic power freq spectral density')
    iret = nf90_put_att(ncid,varid(18),'units','m2 s')
    iret = nf90_put_att(ncid,varid(18),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(18),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(18),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(18),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(18),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(18),'content','TY')
    iret = nf90_put_att(ncid,varid(18),'associates','time frequency')

!edzfth
    iret = nf90_def_var(ncid, 'edzfth_bis',NF90_DOUBLE,(/dimid1(5),dimid1(6),dimid1(1)/),varid(19))

    iret = nf90_put_att(ncid,varid(19),'long_name','ground acoustic directional power spectral density - sum over freq of F_spec for frequency bands ')
    iret = nf90_put_att(ncid,varid(19),'standard_name','ground acoustic directional power spectral density (freq_band)')
    iret = nf90_put_att(ncid,varid(19),'globwave_name','ground acoustic directional power spectral density (freq_band)')
    iret = nf90_put_att(ncid,varid(19),'units','m2 s')
    iret = nf90_put_att(ncid,varid(19),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(19),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(19),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(19),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(19),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(19),'content','TZA')
    iret = nf90_put_att(ncid,varid(19),'associates','direction time frequency')

!edzth
    iret = nf90_def_var(ncid, 'edzth_bis',NF90_DOUBLE,(/dimid1(5),dimid1(1)/),varid(20))

    iret = nf90_put_att(ncid,varid(20),'long_name','ground acoustic power azi spectral density - sum over freq of F_spec')
    iret = nf90_put_att(ncid,varid(20),'standard_name','ground acoustic power azi spectral density')
    iret = nf90_put_att(ncid,varid(20),'globwave_name','ground acoustic power azi spectral density')
    iret = nf90_put_att(ncid,varid(20),'units','m2 s')
    iret = nf90_put_att(ncid,varid(20),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(20),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(20),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(20),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(20),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(20),'content','TA')
    iret = nf90_put_att(ncid,varid(20),'associates','direction time')

!'edzmaz' edz for max of F_spec over AZ
    iret = nf90_def_var(ncid, 'edzmaz_bis',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(21))

    iret = nf90_put_att(ncid,varid(21),'long_name','max over az of F_spec (ground acoustic power spectral density)')
    iret = nf90_put_att(ncid,varid(21),'standard_name','max of F_spec over az')
    iret = nf90_put_att(ncid,varid(21),'globwave_name','max of F_spec over az')
    iret = nf90_put_att(ncid,varid(21),'units','m2 s')
    iret = nf90_put_att(ncid,varid(21),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(21),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(21),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(21),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(21),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(21),'content','TY')
    iret = nf90_put_att(ncid,varid(21),'associates','frequency time')
!'azmax' az corres. to max of F_spec over Az
    iret = nf90_def_var(ncid, 'azmax_bis',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(22))

    iret = nf90_put_att(ncid,varid(22),'long_name','azi of max - max of F_spec over freq')
    iret = nf90_put_att(ncid,varid(22),'standard_name','azi of max')
    iret = nf90_put_att(ncid,varid(22),'globwave_name','azi of max')
    iret = nf90_put_att(ncid,varid(22),'units','m2 s')
    iret = nf90_put_att(ncid,varid(22),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(22),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(22),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(22),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(22),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(22),'content','TY')
    iret = nf90_put_att(ncid,varid(22),'associates','frequency time')
!'edzmf' edz for max of F_spec over f
    iret = nf90_def_var(ncid, 'edzmf_bis',NF90_DOUBLE,(/dimid1(5),dimid1(1)/),varid(23))

    iret = nf90_put_att(ncid,varid(23),'long_name','max over freq of F_spec (ground acoustic power spectral density)')
    iret = nf90_put_att(ncid,varid(23),'standard_name','max of F_spec over freq')
    iret = nf90_put_att(ncid,varid(23),'globwave_name','max of F_spec over freq')
    iret = nf90_put_att(ncid,varid(23),'units','m2 s')
    iret = nf90_put_att(ncid,varid(23),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(23),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(23),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(23),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(23),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(23),'content','TA')
    iret = nf90_put_att(ncid,varid(23),'associates','direction time')
!'fmax' freq corres. to max of F_spec over f
    iret = nf90_def_var(ncid, 'fmax_bis',NF90_DOUBLE,(/dimid1(5),dimid1(1)/),varid(24))

    iret = nf90_put_att(ncid,varid(24),'long_name','freq of max - max of F_spec over freq')
    iret = nf90_put_att(ncid,varid(24),'standard_name','freq of max')
    iret = nf90_put_att(ncid,varid(24),'globwave_name','freq of max')
    iret = nf90_put_att(ncid,varid(24),'units','m2 s')
    iret = nf90_put_att(ncid,varid(24),'scale_factor',1.)
    iret = nf90_put_att(ncid,varid(24),'add_offset',0.)
    iret = nf90_put_att(ncid,varid(24),'valid_min',0.)
    iret = nf90_put_att(ncid,varid(24),'valid_max',1.E20)
    iret = nf90_put_att(ncid,varid(24),'_FillValue',NF90_FILL_DOUBLE)
    iret = nf90_put_att(ncid,varid(24),'content','TA')
    iret = nf90_put_att(ncid,varid(24),'associates','direction time')

!OLD OUTPUTS
    !edzf
!  iret = nf90_def_var(ncid, 'edzf_bis',NF90_DOUBLE,(/dimid1(4),dimid1(1)/),varid(15))

!  iret = nf90_put_att(ncid,varid(15),'long_name',&
!'with alternative Power coeff : ground acoustic power spectral density')
!  iret = nf90_put_att(ncid,varid(15),'standard_name',&
!'with alternative Power coeff : ground acoustic power spectral density')
!  iret = nf90_put_att(ncid,varid(15),'globwave_name',&
!'with alternative Power coeff : ground acoustic power spectral density')
!  iret = nf90_put_att(ncid,varid(15),'units','m2 s')
!  iret = nf90_put_att(ncid,varid(15),'scale_factor',1.)
!  iret = nf90_put_att(ncid,varid(15),'add_offset',0.)
!  iret = nf90_put_att(ncid,varid(15),'valid_min',0.)
!  iret = nf90_put_att(ncid,varid(15),'valid_max',1.E20)
!  iret = nf90_put_att(ncid,varid(15),'_FillValue',NF90_FILL_DOUBLE)
!  iret = nf90_put_att(ncid,varid(15),'content','TY')
!  iret = nf90_put_att(ncid,varid(15),'associates','time frequency')

!edfth
!  iret = nf90_def_var(ncid, 'edfth_bis',NF90_DOUBLE,(/dimid1(5),dimid1(4),dimid1(1)/),varid(16))

!  iret = nf90_put_att(ncid,varid(16),'long_name',&
!'with alternative Power coeff : ground acoustic directional power spectral density')
!  iret = nf90_put_att(ncid,varid(16),'standard_name',&
!'with alternative Power coeff : ground acoustic directional power spectral density')
!  iret = nf90_put_att(ncid,varid(16),'globwave_name',&
!'with alternative Power coeff : ground acoustic directional power spectral density')
!  iret = nf90_put_att(ncid,varid(16),'units','m2 s')
!  iret = nf90_put_att(ncid,varid(16),'scale_factor',1.)
!  iret = nf90_put_att(ncid,varid(16),'add_offset',0.)
!  iret = nf90_put_att(ncid,varid(16),'valid_min',0.)
!  iret = nf90_put_att(ncid,varid(16),'valid_max',1.E20)
!  iret = nf90_put_att(ncid,varid(16),'_FillValue',NF90_FILL_DOUBLE)
!  iret = nf90_put_att(ncid,varid(16),'content','TYA')
!  iret = nf90_put_att(ncid,varid(16),'associates','direction time frequency')

!dist_max
!  iret = nf90_def_var(ncid, 'dist_source_max_bis',NF90_DOUBLE,(/dimid1(5),dimid1(4),dimid1(1)/),varid(17))

!  iret = nf90_put_att(ncid,varid(17),'long_name',&
!'with alternative Power coeff : distance of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(17),'standard_name',&
!'with alternative Power coeff : distance of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(17),'globwave_name',&
!'with alternative Power coeff : distance of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(17),'units','m2 s')
!  iret = nf90_put_att(ncid,varid(17),'scale_factor',1.)
!  iret = nf90_put_att(ncid,varid(17),'add_offset',0.)
!  iret = nf90_put_att(ncid,varid(17),'valid_min',0.)
!  iret = nf90_put_att(ncid,varid(17),'valid_max',1.E20)
!  iret = nf90_put_att(ncid,varid(17),'_FillValue',NF90_FILL_DOUBLE)
!  iret = nf90_put_att(ncid,varid(17),'content','TYA')
!  iret = nf90_put_att(ncid,varid(17),'associates','direction time frequency')

!amp_max for BIS
!  iret = nf90_def_var(ncid, 'source_max_bis',NF90_DOUBLE,(/dimid1(5),dimid1(4),dimid1(1)/),varid(18))

!  iret = nf90_put_att(ncid,varid(18),'long_name',&
!'with alternative Power coeff : amplitude of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(18),'standard_name',&
!'with alternative Power coeff : amplitude of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(18),'globwave_name',&
!'with alternative Power coeff : amplitude of the maximum source (after attenuation) in each direction')
!  iret = nf90_put_att(ncid,varid(18),'units','m2 s')
!  iret = nf90_put_att(ncid,varid(18),'scale_factor',1.)
!  iret = nf90_put_att(ncid,varid(18),'add_offset',0.)
!  iret = nf90_put_att(ncid,varid(18),'valid_min',0.)
!  iret = nf90_put_att(ncid,varid(18),'valid_max',1.E20)
!  iret = nf90_put_att(ncid,varid(18),'_FillValue',NF90_FILL_DOUBLE)
!  iret = nf90_put_att(ncid,varid(18),'content','TYA')
!  iret = nf90_put_att(ncid,varid(18),'associates','direction time frequency')
  END IF

! Add values in NetCDF file

  iret = nf90_enddef(ncid)

  RETURN

END SUBROUTINE W3CRNC
!---------------------------------------------------------------------

!---------------------------------------------------------------------
SUBROUTINE check_err(iret,msg)
  USE NETCDF
  implicit none
  
  integer, INTENT(IN)           :: iret
  CHARACTER(len=*), INTENT(IN)  :: MSG

  IF (iret.NE. NF90_NOERR) THEN
    print *, nf90_strerror(iret)
    print *, msg
    stop
  END IF
  RETURN

END SUBROUTINE check_err
!---------------------------------------------------------------------


! PLAN DU CODE
! Lire inputs .inp
! 
! lire les fichiers qui ne vont pas bouger : Pmap, lat, lon
! 
! Par station : lire les fichiers de sommation...
! Lire le fichier de vent : quand ? (tout lire ou lire jour par jour? A OPTIMISER)
!    Par dt : 
!		- lire le fichier de source
!		- multiplier source par Pmap
!		- boucle sur tous les azimuts :  
!			* boucle sur toutes les cases à sommer, 
!			* calculer attenuation puis 
!			* sommer
!			* ranger la distance du max dans une variable (pour pouvoir situer la source dominante par azimuth)
!			* Pour la distance : stocker min et max dist de la "tache" (+ d'une tache ?)
! 			
! 
! Ranger dans la variable (ou stocker au fur et à mesure ?) Amp à station (dépendant de theta, frequ, temps)
! Sauver dans .nc
! 
!  














