  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---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d
   USE cometdyn
   USE cometemp 
   USE 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





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

   !PRINT*, 'enter METEO routine, nbouc =',nbouc

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

   IF(nbouc==0) THEN

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

! atmospheric pressure
!-----------------------

     pa(:,:)=paref


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



       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
       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



       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)

       ventvent(:,:) = wind_veloc

     ENDIF

! case of unstationary and homogenous wind
! ----------------------------------------
     IF(.NOT.l_meteo_stat .AND. l_meteo_hom) THEN
       IF (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
       ALLOCATE(uvent(nvalmax),vvent(nvalmax))


       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


       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)


       tint1 = tool_datosec(rec(kk:kk+18))
       t1=tchrono(i)
       DO WHILE(t1 < tint1.AND.i<=nvalmax)
         uvent(i)=val11
         vvent(i)=val12

         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


             READ(rec(kk+19:),*) spd,cap

             IF ((spd >= 0.0_rsh).AND.(spd < 100.0_rsh).AND. &
                 (cap >= 0.0_rsh).AND.(cap < 360.0_rsh)) THEN


               val21=-spd*SIN((cap-grid_angle)*pi/180.0_rsh)
               val22=-spd*COS((cap-grid_angle)*pi/180.0_rsh)

               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


                 i=i+1
                 IF(i > nvalmax) GOTO 999
               END DO
               tint1=tint2
               val11=val21
               val12=val22



             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

       ENDIF

       CLOSE(29)
     ENDIF

   END IF ! end test on nbouc

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

   IF(l_meteo_stat) THEN

       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)

       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)




     tosx(:,:)    =dt1*tosu2+(1.0_rlg-dt1)*tosu1
     tosy(:,:)    =dt1*tosv2+(1.0_rlg-dt1)*tosv1

     ventvent(:,:)=dt1*wind2+(1.0_rlg-dt1)*wind1
     


   END IF


   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


   !PRINT*, 'exit METEO routine'

  END SUBROUTINE meteo
