  SUBROUTINE meteo

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE flx_meteo  **
   !&E
   !&E ** Purpose : compute meteorological conditions in case of l_meteo_hom=.true.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : main, step
   !&E
   !&E ** External calls : tool_datosec
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables : tosx, tosy, ventvent, pa, rh, sat, fcc
   !&E                         in/out : uvent, vvent
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-11  (S. Michel, HOCER) Check wind data units and file existence.
   !&E       !  2017-08 (A. Fiandrino, B. Thouvenin. ) introducton of specificities key_lagoon 
   !&E       
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters
   USE comvars2d
   USE cometdyn
   USE cometemp 
   USE_MPI toolmpi, ONLY: MASTER
!#ifdef key_lagoon
!   USE comfluxdo,  ONLY : rain 
!#endif

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   LOGICAL               :: ex
   CHARACTER(len=lchain) :: rec
   CHARACTER(LEN=19)     :: tool_sectodat,date
   INTEGER               :: i,kk,eof
   REAL(kind=rsh)        :: val11,val21,val12,val22,cos1,cos2,spd,cap
   REAL(kind=rsh)        :: wind1,wind2,tosu1,tosv1,tosu2,tosv2
   REAL(kind=rlg)        :: dt1,tint1,tint2,tool_datosec,t1

#ifdef key_lagoon
   INTEGER j
   REAL(kind=rsh)        :: patm_lu1,tempair_lu1,relh_lu1,nebu_lu1
   REAL(kind=rsh)        :: patm_lu2,tempair_lu2,relh_lu2,nebu_lu2
   REAL(kind=rsh)        :: val1,val2,val3,val1b,val2b,val3b
#endif

#ifdef key_lagoon_spec
! le 28/05/2013 Abattement du vent dans les zones conch de la lagune de Thau
   REAL(kind=rsh)        :: abat1,abat2
   REAL(kind=rsh)        :: val11_abat,val21_abat,val12_abat,val22_abat
   REAL(kind=rsh)        :: wind1_abat,wind2_abat,tosu1_abat,tosv1_abat
   REAL(kind=rsh)        :: tosu2_abat,tosv2_abat
#endif

   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter METEO routine, nbouc =',nbouc

   !  initialisation, reading data from the input file
   !--------------------------------------------------

   IF(nbouc==0) THEN

     ! save into simu.log
     !-------------------
     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) '**************** METEO.F90 ***********************'
     WRITE(iscreenlog,*) '**************************************************'
     ENDIF_MPI

     ! atmospheric pressure
     !-----------------------
#if ! defined key_equeta_temp
     pa(:,:)=paref
#endif

     ! case of stationary and homogenous wind
     !---------------------------------------
     IF(l_meteo_stat) THEN
       IF_MPI (MASTER) THEN
       WRITE(iscreenlog,*) ' '

#ifdef key_lagoon_spec
! le 28/05/2013 Abattement du vent dans les zones conch de la lagune de Thau
       IF (wind_dir >= 11.25_rsh .AND. wind_dir <123.75_rsh) abat1=0.38_rsh
       IF (wind_dir >= 123.75_rsh .AND. wind_dir < 213.75_rsh) abat1=0.27_rsh
       IF (wind_dir >= 213.75_rsh .AND. wind_dir < 281.25_rsh) abat1=0.38_rsh
       IF (wind_dir >= 281.25_rsh .OR. wind_dir < 11.25_rsh) abat1=0.59_rsh
#endif

       wind_dir = MODULO(wind_dir-grid_angle,360.0_rsh)    

       WRITE(iscreenlog,*) 'Stationary and homogeneous wind'
       WRITE(iscreenlog,*) 'velocity (in m.s-1): wind_veloc = ',wind_veloc
       WRITE(iscreenlog,*) 'direction (0°: from North): wind_dir = ',wind_dir
       ENDIF_MPI
       IF ((wind_veloc<0.0_rsh).OR.(wind_veloc>100.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In METEO.F90, the constant wind velocity'
         WRITE(ierrorlog,*) 'is not expressed in m.s-1 :'
         WRITE(ierrorlog,*) 'wind_veloc = ',wind_veloc
         WRITE(ierrorlog,*) 'Correct the "paracom.txt" parameters file.'
         WRITE(ierrorlog,*) 'Simulation stopped.'
         STOP
       ENDIF
       IF ((wind_dir<0.0_rsh).OR.(wind_dir>=360.0_rsh)) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In METEO.F90, the constant wind direction'
         WRITE(ierrorlog,*) 'is not expressed in degrees from the northerly direction:'
         WRITE(ierrorlog,*) 'wind_dir = ',wind_dir
         WRITE(ierrorlog,*) 'Correct the "paracom.txt" parameters file.'
         WRITE(ierrorlog,*) 'Simulation stopped.'
         STOP
       ENDIF

#ifdef key_lagoon
#ifdef key_lagoon_spec
       tosx(:,:) = -(1-xmask_hx(:,:) + xmask_hx(:,:)*abat1) &
                   *wind_veloc*wind_veloc*SIN(wind_dir*pi/180.0_rsh)
       tosy(:,:) = -(1-xmask_hy(:,:) + xmask_hy(:,:)*abat1) &
                   *wind_veloc*wind_veloc*COS(wind_dir*pi/180.0_rsh)
       ventvent(:,:) = (1-xmask_h0(:,:) + xmask_h0(:,:)*abat1)*wind_veloc
#else
       tosx(:,:) = -wind_veloc*wind_veloc*SIN(wind_dir*pi/180.0_rsh)
       tosy(:,:) = -wind_veloc*wind_veloc*COS(wind_dir*pi/180.0_rsh)     
       ventvent(:,:) = wind_veloc
#endif    
#else  
#if ! defined key_palm_ww3
       tosx(:,:) = -wind_veloc*wind_veloc*SIN((wind_dir-grid_angle)*pi/180.0_rsh)
       tosy(:,:) = -wind_veloc*wind_veloc*COS((wind_dir-grid_angle)*pi/180.0_rsh)
#endif
       ventvent(:,:) = wind_veloc
#endif
     ENDIF

     ! case of unstationary and homogenous wind
     ! ----------------------------------------
     IF(.NOT.l_meteo_stat .AND. l_meteo_hom) THEN
       IF_MPI (MASTER) THEN
       WRITE(iscreenlog,*) ' '
       WRITE(iscreenlog,*) 'Unstationary and homogeneous wind'
       WRITE(iscreenlog,*) 'Time-series of wind speed and direction are read from the file:'
       WRITE(iscreenlog,*) TRIM(file_meteo)
       ENDIF_MPI
       ALLOCATE(uvent(nvalmax),vvent(nvalmax))
#ifdef key_lagoon
       ALLOCATE(patm(nvalmax),tempair(nvalmax),relh(nvalmax),nebu(nvalmax))
#ifdef key_lagoon_spec
       ALLOCATE(uvent_abat(nvalmax),vvent_abat(nvalmax))
#endif
#endif

       INQUIRE(file=file_meteo,exist=ex)
       IF (.NOT.ex) THEN
         WRITE(ierrorlog,*) ' '
         WRITE(ierrorlog,*) 'ERROR: In METEO.F90, wind input time-series are needed '
         WRITE(ierrorlog,*) 'to compute unstationary wind stress.'
         WRITE(ierrorlog,*) 'Check the path and name of the meteorological file:'
         WRITE(ierrorlog,*) TRIM(file_meteo)
         WRITE(ierrorlog,*) 'or set l_meteo_stat=TRUE to use a constant wind field.'
         WRITE(ierrorlog,*) 'Simulation stopped.'
         STOP
       ENDIF

       OPEN(29,file = file_meteo,form='formatted')
       i=1

#ifdef key_lagoon
       READ(29,'(a)',END=999,err=999) rec
       kk=INDEX(rec,'/')-2
       IF (l_cmeteocte) THEN
         READ(rec(kk+19:),*) spd,cap,patm_lu1,tempair_lu1
         relh_lu1 = rhcte ! valeur en %
         nebu_lu1 = nebcte
         relh_lu2 = rhcte
         nebu_lu2 = nebcte      
       ELSE
         READ(rec(kk+19:),*) spd,cap,patm_lu1,tempair_lu1,relh_lu1,nebu_lu1
! conversion des valeurs lues comprises entre 0 et 1 
         relh_lu1=relh_lu1/100_rsh    
         nebu_lu1=nebu_lu1/8.0_rsh 
       ENDIF          

#ifdef key_lagoon_spec
! le 28/05/2013 Abattement du vent dans les zones conch de la lagune de Thau
       IF (cap >= 11.25_rsh .AND. cap < 123.75_rsh) abat1=0.38_rsh
       IF (cap >= 123.75_rsh .AND. cap < 213.75_rsh) abat1=0.27_rsh
       IF (cap >= 213.75_rsh .AND. cap < 281.25_rsh) abat1=0.38_rsh
       IF (cap >= 281.25_rsh .OR. cap <11.25_rsh) abat1=0.59_rsh
#endif

       cap=MODULO(cap-grid_angle,360.0_rsh)
       
       val11=-spd*SIN(cap*pi/180.0_rsh)
       val12=-spd*COS(cap*pi/180.0_rsh)

#ifdef key_lagoon_spec
! le 28/05/2013 Abattement du vent dans les zones conch de la lagune de Thau
       val11_abat=-abat1*spd*SIN(cap*pi/180.0_rsh)
       val12_abat=-abat1*spd*COS(cap*pi/180.0_rsh) 
#endif

#else
       READ(29,*)
       READ(29,*)
       READ(29,*)
       READ(29,*)
       READ(29,'(a)',END=999,err=999) rec
       kk=INDEX(rec,'/')-2
       READ(rec(kk+19:),*) spd,cap
       val11=-spd*SIN((cap-grid_angle)*pi/180.0_rsh)
       val12=-spd*COS((cap-grid_angle)*pi/180.0_rsh)
#endif

       tint1 = tool_datosec(rec(kk:kk+18))
       t1=tchrono(i)
       DO WHILE(t1 < tint1.AND.i<=nvalmax)
         uvent(i)=val11
         vvent(i)=val12
#ifdef key_lagoon
         patm(i)=patm_lu1
         tempair(i)=tempair_lu1
         relh(i)=relh_lu1
         nebu(i)=nebu_lu1
#endif
         i=i+1
         IF(i<=nvalmax) t1=tchrono(i)
       END DO

       ! If the simulation starts before the first data record,
       ! the wind velocities are set to the first available data.
       IF(i > 1) THEN
         date=tool_sectodat(tint1)
         WRITE(iwarnlog,*) ' '
         WRITE(iwarnlog,*) 'WARNING: In METEO.F90, the wind time series'
         WRITE(iwarnlog,*) 'is not available from the beginning of the simulation.'
         WRITE(iwarnlog,*) 'The first values from the file: ',TRIM(file_meteo)
         WRITE(iwarnlog,*) 'corresponding to the date : ',date
         WRITE(iwarnlog,*) 'are repeated for the first ',i,' input time-steps.'
       ENDIF

       ! time interpolation of wind data
       IF(i < nvalmax) THEN
         eof=0
         DO WHILE(eof==0) 
           READ(29,'(a)',iostat=eof) rec
           IF(eof==0) THEN
             kk=INDEX(rec,'/')-2

#ifdef key_lagoon
             IF (l_cmeteocte) THEN
               READ(rec(kk+19:),*) spd,cap,patm_lu2,tempair_lu2
             ELSE
               READ(rec(kk+19:),*) spd,cap,patm_lu2,tempair_lu2,relh_lu2,nebu_lu2
               relh_lu2=relh_lu2/100_rsh
               nebu_lu2=nebu_lu2/8.0_rsh
             ENDIF     
#else
             READ(rec(kk+19:),*) spd,cap
#endif
             IF ((spd >= 0.0_rsh).AND.(spd < 100.0_rsh).AND. &
                 (cap >= 0.0_rsh).AND.(cap < 360.0_rsh)) THEN

#ifdef key_lagoon
#ifdef key_lagoon_spec
! le 28/05/2013 Abattement du vent dans les zones conch de la lagune de Thau
               IF (cap >= 11.25_rsh .AND. cap < 123.75_rsh) abat2=0.38_rsh
               IF (cap >= 123.75_rsh .AND. cap < 213.75_rsh) abat2=0.27_rsh
               IF (cap >= 213.75_rsh .AND. cap < 281.25_rsh) abat2=0.38_rsh
               IF (cap >= 281.25_rsh .OR. cap <11.25_rsh) abat2=0.59_rsh
#endif
               cap=MODULO(cap-grid_angle,360.0_rsh)
               val21=-spd*SIN(cap*pi/180.0_rsh)
               val22=-spd*COS(cap*pi/180.0_rsh)

#ifdef key_lagoon_spec
               val21_abat=-abat2*spd*SIN(cap*pi/180.0_rsh)
               val22_abat=-abat2*spd*COS(cap*pi/180.0_rsh) 
#endif      

#else
               val21=-spd*SIN((cap-grid_angle)*pi/180.0_rsh)
               val22=-spd*COS((cap-grid_angle)*pi/180.0_rsh)
#endif
               tint2 = tool_datosec(rec(kk:kk+18))
               DO WHILE ((tchrono(i)-tint1)*(tchrono(i)-tint2) <= 0.0_rlg)
                 dt1=(tchrono(i)-tint1)/(tint2-tint1)
                 uvent(i)=val11*(1.0_rlg-dt1)+val21*dt1
                 vvent(i)=val12*(1.0_rlg-dt1)+val22*dt1

#ifdef key_lagoon
                 patm(i)=patm_lu1*(1.d0-dt1)+patm_lu2*dt1
                 tempair(i)=tempair_lu1*(1.d0-dt1)+tempair_lu2*dt1
                 relh(i)=relh_lu1*(1.d0-dt1)+relh_lu2*dt1            
                 nebu(i)=nebu_lu1*(1.d0-dt1)+nebu_lu2*dt1
#ifdef key_lagoon_spec
                 uvent_abat(i)=val11_abat*(1.0_rlg-dt1)+val21_abat*dt1
                 vvent_abat(i)=val12_abat*(1.0_rlg-dt1)+val22_abat*dt1
#endif

#endif
                 i=i+1
                 IF(i > nvalmax) GOTO 999
               END DO
               tint1=tint2
               val11=val21
               val12=val22
#ifdef key_lagoon_spec
               val11_abat=val21_abat
               val12_abat=val22_abat
#endif

#ifdef key_lagoon
               patm_lu1=patm_lu2
               tempair_lu1=tempair_lu2
               relh_lu1=relh_lu2
               nebu_lu1=nebu_lu2
#endif
             ELSE
               WRITE(ierrorlog,*) ' '
               WRITE(ierrorlog,*) 'ERROR: In METEO.F90, the unstationary wind data must be'
               WRITE(ierrorlog,*) 'velocity expressed in m.s-1 and  '
               WRITE(ierrorlog,*) 'direction expressed in degrees from the northerly direction.'
               WRITE(ierrorlog,*) 'Check the input file containing the wind time-series: '
               WRITE(ierrorlog,*) TRIM(file_meteo)
               WRITE(ierrorlog,*) 'record number = ',i
               WRITE(ierrorlog,*) 'speed = ',spd,', direction = ',cap
               WRITE(ierrorlog,*) 'Simulation stopped.'
               STOP
             END IF
           END IF
         END DO
       END IF

999    CONTINUE

       ! If the simulation ends after the last data record,
       ! set the wind velocities to the last available data.
       IF(i < nvalmax) THEN
         date=tool_sectodat(tint2)
         WRITE(iwarnlog,*) ' '
         WRITE(iwarnlog,*) 'WARNING: In METEO.F90, the wind time series'
         WRITE(iwarnlog,*) 'is not available to the end of the simulation.'
         WRITE(iwarnlog,*) 'The last values from the file: ',TRIM(file_meteo)
         WRITE(iwarnlog,*) 'corresponding to the date : ',date
         WRITE(iwarnlog,*) 'are repeated for the last ',nvalmax-i+1,' input time-steps.'
         uvent(i:nvalmax)=val21
         vvent(i:nvalmax)=val22
#ifdef key_lagoon_spec
         uvent_abat(i:nvalmax)=val21_abat
         vvent_abat(i:nvalmax)=val22_abat
#endif
       ENDIF

       CLOSE(29)
     ENDIF

   END IF ! end test on nbouc

   !**********************
   ! ordinary iteration
   !**********************

   IF(l_meteo_stat) THEN
#if ! defined key_palm_ww3
       tosx(:,:) = -wind_veloc*wind_veloc*SIN((wind_dir-grid_angle)*pi/180.0_rsh)
       tosy(:,:) = -wind_veloc*wind_veloc*COS((wind_dir-grid_angle)*pi/180.0_rsh)
#endif
       ventvent(:,:) = wind_veloc
   ELSE IF(.NOT.l_meteo_stat.AND.l_meteo_hom) THEN
     dt1=(t-tchrono(indchrono))/(tchrono(indchrono+1)-tchrono(indchrono))
     wind1=SQRT(uvent(indchrono  )**2+vvent(indchrono  )**2)
     wind2=SQRT(uvent(indchrono+1)**2+vvent(indchrono+1)**2)
     tosu1=wind1*uvent(indchrono  )
     tosv1=wind1*vvent(indchrono  )
     tosu2=wind2*uvent(indchrono+1)
     tosv2=wind2*vvent(indchrono+1)

#ifdef key_lagoon_spec
     wind1_abat=SQRT(uvent_abat(indchrono  )**2+vvent_abat(indchrono  )**2)
     wind2_abat=SQRT(uvent_abat(indchrono+1)**2+vvent_abat(indchrono+1)**2)
     tosu1_abat=wind1_abat*uvent_abat(indchrono  )
     tosv1_abat=wind1_abat*vvent_abat(indchrono  )
     tosu2_abat=wind2_abat*uvent_abat(indchrono+1)
     tosv2_abat=wind2_abat*vvent_abat(indchrono+1)
#endif

#if ! defined key_palm_ww3
     tosx(:,:)    =dt1*tosu2+(1.0_rlg-dt1)*tosu1
     tosy(:,:)    =dt1*tosv2+(1.0_rlg-dt1)*tosv1
#endif
     ventvent(:,:)=dt1*wind2+(1.0_rlg-dt1)*wind1
     
#ifdef key_lagoon
#ifdef key_lagoon_spec
     DO j=ljmin,ljmax
        tosx(limin:limax,j)    =dt1*((1._rsh-xmask_hx(limin:limax,j))*tosu2+xmask_hx(limin:limax,j)*tosu2_abat) &
                    +(1.0_rlg-dt1)*((1._rsh-xmask_hx(limin:limax,j))*tosu1+xmask_hx(limin:limax,j)*tosu1_abat)
        tosy(limin:limax,j)    =dt1*((1._rsh-xmask_hy(limin:limax,j))*tosv2+xmask_hy(limin:limax,j)*tosv2_abat) &
                    +(1.0_rlg-dt1)*((1._rsh-xmask_hy(limin:limax,j))*tosv1+xmask_hy(limin:limax,j)*tosv1_abat)
        ventvent(limin:limax,j) =dt1*((1._rsh-xmask_h0(limin:limax,j))*wind2+xmask_h0(limin:limax,j)*wind2_abat) &
                   +(1_rlg-dt1)*((1._rsh-xmask_h0(limin:limax,j))*wind1+xmask_h0(limin:limax,j)*wind1_abat)
     END DO
#endif
     pa(:,:)=dt1*patm(indchrono+1)+(1.d0-dt1)*patm(indchrono)
     sat(:,:)=dt1*tempair(indchrono+1)+(1.d0-dt1)*tempair(indchrono)
     rh(:,:)=dt1*relh(indchrono+1)+(1.d0-dt1)*relh(indchrono)
     fcc(:,:)=dt1*nebu(indchrono+1)+(1.d0-dt1)*nebu(indchrono)
#endif

   END IF

#if ! defined key_lagoon
   cos1=COS(2.0_rsh*pi/365.25_rsh*REAL(jjulien))
   cos2=COS(2.0_rsh*pi/365.25_rsh*(REAL(jjulien)-220.0_rsh))
   IF(l_meteo_hom) fcc(:,:)=(65.0_rsh+15.0_rsh*cos1)/100.0_rsh  ! 0 < fcc < 1
   IF(l_meteo_hom) sat(:,:)=13.0_rsh+4.0_rsh*cos2
   IF(l_meteo_hom) rh(:,:) =80.0_rsh/100.0_rsh                  ! 0 < rh < 1
#endif

   PRINT_DBG*, 'exit METEO routine'

  END SUBROUTINE meteo
