  SUBROUTINE sflx_surf(temp,sal,xe3)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sflx_surf  ***
   !&E
   !&E ** Purpose : call the routines computing exchanges at the air-sea interface
   !&E              (except for momentum exchanges)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step3dxy,step3dyx,step2dxy,step2dyx
   !&E
   !&E
   !&E ** Used ij-arrays : temp,sal,xe3,h0,surf,tclim
   !&E                     rad,irlos,senslos,latlos,dirlos,dlatlos,dsenslos,hfrelax,xe3,
   !&E                     extinction,phitemp,phisal,phiriver,phieau,rain,evap,hfrelax,
   !&E                     obc_ts,obc_ss,obc_tn,obc_sn,obc_tw,obc_sw,obc_te,obc_se
   !&E  
   !&E ** Modified variables : hfrelax,phitemp,phieau,
   !&E                         latlos,dlatlos,evap (sflx_turb),
   !&E                         irlos,dirlos (sflx_ir),
   !&E                         senslos,dsenslos (sflx_turb),
   !&E                         rain (rain)
   !&E
   !&E ** History :
   !&E       !                        Original code
   !&E       !  2008-06  (V. Garnier) introduction heat flux relative
   !&E                   to river discharges and rain-evap flux
   !&E       !  2008-06  (P. Garreau) extinction coefficient relative to depth
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E       !  2009-01  (M. Jouan, V. Garnier) initialisation of phitemp
   !&E                                          for constant radiation (l_sfls_radcst=T)
   !&E       !  2009-09  (B. Thouvenin) addition of atmospheric flux for substances
   !&E       !  2009-11  (S. Michel, HOCER) Check needed met. variables are read. 
   !&E       !  2010-01  (I. Pairaud, P. Garreau, R. Ramel) correction of tempfle estimate
   !&E       !  2010-03  (V. Garnier) correction of phiriver to allow several river
   !&E       !                        discharges at the same location
   !&E       !  2010-10  (B. Thouvenin) Distribution of the river discharge throughout a thickness
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2014-05 (B.Thouvenin) substance atmospheric inputs variable in space and time + zerogradient and OBC relaxation 
   !&E       !  2014-08 (P. Garreau) Introduction of 2 wave length in the formulation of heat content penetration. IMPORTANT modification.
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used

   USE parameters, ONLY : imin,imax,kmax,jmin,jmax,rsh,rlg,pi,num_testcase,obc_width, &
                          limin,liminp1,limax,ljminp1,ljmax,ljmin,valmanq,liminm2,    &
                          liminm1,limaxp2,ljminm1,ljmaxp2,limaxp1,ljmaxp1,ljminm2,l_testcase
   USE comvars2d,  ONLY : h0fond,halfdt,jjulien,iheure,iminu,isec,h0,surf,  &
                          l_modele2d,jb,jh,ig,id,jhemisph,ierrorlog,fwetp,h0_g
   USE comvars3d,  ONLY : rhoref,dsigu
   USE comfluxdo,  ONLY : nbfl,tmfle,iriv,jriv,phieau,phiriver,evap,rain,amtfle,l_tempriv
   USE cometemp,   ONLY : l_rapflux,l_sflx_solarcst,l_sflx_radlossbot,   &
                          phitemp,phisal,tclim,hfrelax,coext,chp,sflx_solarcst,&
                          sflx_rad_maxdepth
   USE sflxatm,    ONLY : dlatlos, rad, irlos, senslos, latlos
   USE sflxrad,    ONLY : sflx_ir, sflx_solar
   USE sflxturb,   ONLY : sflx_turb, lnu
   USE obccombine, ONLY : l_obc_south, l_obc_north, l_obc_west, l_obc_east
   USE obcogcm,    ONLY : l_obc_ogcm_relax, obc_coefrel,  &
                          obc_ts,obc_ss,obc_tn,obc_sn,obc_tw, obc_sw, obc_te, obc_se
#if defined key_substance 
   USE comsubstance,ONLY : nv_adv,phicon,phicon_atm,sub_flx_atm,cv_rain, &
                           sflx_sub_atm_depth,l_subflxatm,l_subflxatm_xyt,cv_wat, &
                           obc_cv_wat_n,obc_cv_wat_s,obc_cv_wat_w,obc_cv_wat_e
   USE subreadflxatm,ONLY : sub_flx_atm_xyt
#if defined key_biolo
   USE comvars2d,   ONLY : t
   USE comsubstance,ONLY : extinction,tdebsubs
   USE bioloinit,   ONLY : l_bioretro_extinct
#endif
#if defined key_zerograd
   USE comsubstance,ONLY : obc_cv_name
#endif
#endif
   USE toolgeom,    ONLY : f_zw,f_dzu
#if defined key_daily_climato_kpar
   USE comvars2d,   ONLY : nbouc,kpar_sat,jjulien_last,iscreen
   USE ionc4,       ONLY : ionc4_openr, ionc4_read_subxyt, ionc4_close
   USE_MPI toolmpi, ONLY : MASTER
#endif

   IMPLICIT NONE

   !! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(IN) :: xe3
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(IN) :: sal
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(IN) :: temp

   !! * Local declarations
   INTEGER         :: i,j,k,iv,ncerror
   REAL(kind=rsh)  :: d3,rad_top,rad_bot,tempfle,dz,immer,dz3,immer_top,immer_bot
   REAL(kind=rlg)  :: tjour,phi
   REAL(kind=rsh)  :: ms,eta1,eta2

#if defined key_daily_climato_kpar
   CHARACTER(len=200) :: climato_file
   REAL(kind=rsh), DIMENSION(limin:limax,ljmin:ljmax) :: tmp_2d_rsh
#endif

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

!eau tres claire
!  ms=0.40
!  eta1=0.20
!  eta2=0.025
!eau  claire
   ms=0.62_rsh
   eta1=0.66_rsh
   eta2=0.05_rsh
!eau turbide
!  ms=0.78_rsh
!  eta1=0.66_rsh
!  eta2=0.125_rsh
! nemo
!   ms =  0.58_rsh
!   eta1 = 1./0.35_rsh
!   eta2 = 1./23_rsh


   IF (kmax/=1.AND.l_modele2d) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'ERROR: routine SFLXSURF.F90'
     PRINT*,'l_modele2d = .true. and kmax differs from 1'
     PRINT*,'The calculation of "phitemp" is uncorrect in these conditions:'
     PRINT*,'l_modele2d = .true. , kmax /= 1'
     PRINT*,'Simulation stopped'
     STOP
   END IF

#if defined key_daily_climato_kpar
   climato_file = '../../inputs/KPAR/clim_kpar_daily_BREST50m.nc'
   IF (nbouc == 1) THEN
     jjulien_last = 0
     kpar_sat(:,:)=-valmanq
   END IF

   IF (jjulien_last /= jjulien) THEN
     jjulien_last = jjulien

     IF_MPI (MASTER) THEN
       WRITE(iscreen,*) 'Update of daily kpar sat',jjulien_last
     ENDIF_MPI

     CALL ionc4_openr(trim(climato_file))
     CALL ionc4_read_subxyt(trim(climato_file),'kpar',tmp_2d_rsh,limin,limax,ljmin,ljmax,jjulien_last,0,0)
     kpar_sat(limin:limax,ljmin:ljmax) = max(tmp_2d_rsh(limin:limax,ljmin:ljmax),0.1_rsh)
     CALL ionc4_close(trim(climato_file))
   END IF
#endif

!$OMP PARALLEL DEFAULT(NONE) SHARED(l_rapflux,l_sflx_solarcst           &
!$OMP& ,sflx_rad_maxdepth,temp,hfrelax,dlatlos                          &
#if defined key_MPI_2D || defined key_substance
!$OMP& ,limin                                                           &
#endif
!$OMP& ,surf,liminp1,limax,ljmin,ljmax,halfdt                           &
!$OMP& ,tjour,nbfl,tempfle,iriv,jriv,phitemp,phieau,phiriver            &
!$OMP& ,l_obc_ogcm_relax,l_obc_south,l_obc_north,l_obc_west,l_obc_east  &
!$OMP& ,jjulien,iheure,iminu,isec,l_tempriv,tmfle,amtfle,tclim,h0,h0_g  &
!$OMP& ,xe3,rad,irlos,senslos,latlos,obc_ts,phisal,obc_ss,sal           &
!$OMP& ,obc_tn,obc_sn,obc_tw,obc_sw,obc_te,obc_se,rain,evap             &
!$OMP& ,jb,jh,ig,id,jhemisph,fwetp,ms,eta1,eta2                         &
#if defined key_siggen || defined key_gencoord
#else
!$OMP& ,dsigu                                                           &
#endif
#ifdef key_substance
!$OMP & ,phicon,phicon_atm,sub_flx_atm,cv_rain,nv_adv,l_subflxatm       &
!$OMP & ,sflx_sub_atm_depth,cv_wat,l_subflxatm_xyt,sub_flx_atm_xyt      &
!$OMP & ,obc_cv_wat_n,obc_cv_wat_s,obc_cv_wat_w,obc_cv_wat_e            &
#ifdef key_biolo
!$OMP & ,t,tdebsubs,extinction,l_bioretro_extinct,coext                 &
#endif
#ifdef key_zerograd
!$OMP & ,obc_cv_name                                                    &
#endif
#endif
#ifdef key_daily_climato_kpar
!$OMP & ,kpar_sat &
#endif
!$OMP& ,ljminp1,h0fond,l_sflx_radlossbot,chp,rhoref,obc_coefrel,sflx_solarcst) &
#ifdef key_agrif
!$OMP& SHARED(imin,imax,jmin,jmax)                                      &
#endif
!$OMP& private(k,d3,rad_top,rad_bot,i,j,phi,immer,immer_top,iv,immer_bot,dz,dz3,ncerror)

!!!! pour agrif ajouter !$OMP& SHARED(l_testcase,num_testcase) ?

   IF (.NOT.l_sflx_solarcst .OR. (l_testcase .AND.  &
         (num_testcase==36 .OR. num_testcase==37    &
#ifdef key_biolo
                           .OR. num_testcase==40    &
#endif         
           ))) THEN

   ! Applies heat fluxes along the vertical
   !---------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
     DO j = ljminp1,ljmax
       DO i = liminp1,limax
         d3=h0(i,j)+xe3(i,j)
         IF (d3 > h0fond) THEN

           k = kmax
#if defined key_substance && defined key_biolo
           IF (t .GE. tdebsubs .AND. l_bioretro_extinct) coext=MAX(0.154_rsh,extinction(k,i,j))
#endif
           dz = f_dzu(h0(i,j),xe3(i,j),k,i,j)
           immer = dz
           rad_top=rad(i,j)
           !rad_bot=rad_top*EXP(-coext*dz)
#if defined key_daily_climato_kpar
           rad_bot=rad_top*EXP(-kpar_sat(i,j)*dz)
#else
           rad_bot=rad_top*(ms*EXP(-eta1*immer)+(1.0_rsh-ms)*EXP(-eta2*immer))
#endif
           phi=halfdt*(rad_top-rad_bot+irlos(i,j)+senslos(i,j)+latlos(i,j))/(chp*rhoref)
           phitemp(k,i,j) = phi*fwetp(i,j)

           IF(l_rapflux) phitemp(k,i,j)=phitemp(k,i,j)+halfdt*hfrelax(i,j)/(chp*rhoref)*fwetp(i,j)

           DO k = kmax-1,2,-1
#if defined key_substance && defined key_biolo
             IF (t .GE. tdebsubs .AND. l_bioretro_extinct) coext=MAX(0.154_rsh,extinction(k,i,j))
#endif
             ! immersion from the sea surface. The water thickness is xe - z (z being negative)
             immer = - (f_zw(h0(i,j),xe3(i,j),k,i,j)-xe3(i,j))
             dz = f_dzu(h0(i,j),xe3(i,j),k,i,j)
             IF (immer <= sflx_rad_maxdepth) THEN
               rad_top=rad_bot
!               rad_bot=rad_top*EXP(-coext*dz)
#if defined key_daily_climato_kpar
               rad_bot=rad_top*EXP(-kpar_sat(i,j)*dz)
#else
               rad_bot=rad(i,j)*(ms*EXP(-eta1*immer)+(1.0_rsh-ms)*EXP(-eta2*immer))
#endif
               phi=halfdt*(rad_top-rad_bot)/(chp*rhoref)
               phitemp(k,i,j) = phi*fwetp(i,j)
             ELSE
               phitemp(k,i,j)=0.0_rlg
             END IF
           END DO

           k = 1
#if defined key_substance && defined key_biolo
           IF (t .GE. tdebsubs .AND. l_bioretro_extinct) coext=MAX(0.154_rsh,extinction(k,i,j))
#endif
           ! immersion from the sea surface. The water thickness is xe - z (z being negative)
           immer= - (f_zw(h0(i,j),xe3(i,j),k,i,j)-xe3(i,j))
           dz = f_dzu(h0(i,j),xe3(i,j),k,i,j)
           IF (immer <= sflx_rad_maxdepth) THEN
             rad_top=rad_bot
             IF (.NOT.l_sflx_radlossbot) THEN
               phi=halfdt*rad_top/(chp*rhoref)
             ELSE
               !rad_bot=rad_top*EXP(-coext*dz)
#if defined key_daily_climato_kpar
               rad_bot=rad_top*EXP(-kpar_sat(i,j)*dz)
#else
               rad_bot=rad(i,j)*(ms*EXP(-eta1*immer)+(1.0_rsh-ms)*EXP(-eta2*immer))
#endif
               phi=halfdt*(rad_top-rad_bot)/(chp*rhoref)
             ENDIF
             phitemp(k,i,j) = phi*fwetp(i,j)
           ELSE
             phitemp(k,i,j)=0.0_rlg
           ENDIF

         END IF
       END DO
     END DO
!$OMP END DO


   !  water flux due to evaporation precipitation balance
   !  array is in m3/s and evap - rain fields are in mm/s
   !  (1kg/m2 = 10-3m=1mm)
   !  heat flux relative to evaporation precipitation balance
   !-----------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
     DO j=ljminp1,ljmax
       DO i = liminp1,limax
         IF(h0(i,j)+xe3(i,j) > h0fond) THEN

     ! water flux
           phieau(kmax,i,j) = 0.001_rsh*(rain(i,j)-evap(i,j))*surf(i,j)
           IF(l_rapflux) phieau(kmax,i,j) = phieau(kmax,i,j)-                        &
                   0.001_rsh*dlatlos(i,j)/lnu*(temp(kmax,i,j)-tclim(i,j))*surf(i,j)

     ! heat flux
           phitemp(kmax,i,j) = phitemp(kmax,i,j) + halfdt*0.001_rlg*rain(i,j)*temp(kmax,i,j)*fwetp(i,j)

         END IF
       END DO
     END DO
!$OMP END DO


   ! apport de chaleur par les fleuves
   !-----------------------------------

!$OMP SINGLE
     tjour=jjulien+iheure/24.0_rlg+iminu/1440.0_rlg+isec/86400.0_rlg
     DO k = 1,nbfl
       IF (l_tempriv(k)) THEN
         tempfle=tmfle(k)
       ELSE
         tempfle=tmfle(k)+amtfle(k)*REAL(COS(2.0_rlg*pi/365.256363_rlg*(tjour-220.0_rlg)),rsh)*jhemisph
       ENDIF
       i=iriv(k)
       j=jriv(k)
       IF_MPI (ljmin<=j .AND. j<=ljmax) THEN
       IF_MPI (limin<=i .AND. i<=limax) THEN
         phitemp(:,i,j) = phitemp(:,i,j)+halfdt*phiriver(k,:)*tempfle/surf(i,j)
         phieau(1:kmax,i,j)=phieau(1:kmax,i,j) + phiriver(k,1:kmax)
       ENDIF_MPI
       ENDIF_MPI
     END DO
!$OMP END SINGLE


   ! relaxation to T and S and substances along open boundaries
   !--------------------------------------------
     IF (l_obc_ogcm_relax) THEN
!toseelater has not been tested yet with MPI
       IF (l_obc_south) THEN
!$OMP DO SCHEDULE(STATIC)
         DO j = jmin+2,jmin+obc_width
           IF_MPI (ljmin<=j .AND. j<=ljmax) THEN
             DO i = liminp1,limax
               IF ( h0_g(i,jb(i)) /= -valmanq) THEN
#if defined key_siggen || defined key_gencoord
                 DO k=1,kmax
                   dz3=f_dzu(h0(i,j),xe3(i,j),k,i,j)
                   phitemp(k,i,j)=phitemp(k,i,j)-halfdt*obc_coefrel*(temp(k,i,j)-obc_ts(k,i,j))*dz3*fwetp(i,j)
                   phisal(k,i,j) =    0.0_rlg -halfdt*obc_coefrel*(sal(k,i,j)-obc_ss(k,i,j))*dz3*fwetp(i,j)
#ifdef key_substance
#ifdef key_zerograd
                   DO iv=1,nv_adv
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                       phicon(iv,k,i,j) = phicon(iv,k,i,j) -halfdt*obc_coefrel*(cv_wat(iv,k,i,j)-obc_cv_wat_s(iv,k,i,j))*dz3*fwetp(i,j)
                    ENDIF
                  ENDDO
#else
                   phicon(:,k,i,j) = phicon(:,k,i,j) -halfdt*obc_coefrel*(cv_wat(:,k,i,j)-obc_cv_wat_s(:,k,i,j))*dz3*fwetp(i,j)
#endif
#endif
                 END DO
#else
                 phitemp(:,i,j)=phitemp(:,i,j)-halfdt*obc_coefrel*(temp(:,i,j)-obc_ts(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                 phisal(:,i,j) =    0.0_rlg -halfdt*obc_coefrel*(sal(:,i,j)-obc_ss(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#ifdef key_substance
                   DO iv=1,nv_adv
#ifdef key_zerograd
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                       phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_s(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                    ENDIF
#else
                    phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_s(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#endif
                 ENDDO
#endif
#endif
               END IF
             END DO
           ENDIF_MPI
         END DO
!$OMP END DO
       END IF

       IF (l_obc_north) THEN
!$OMP DO SCHEDULE(STATIC)
         DO j = jmax-obc_width+1,jmax
           IF_MPI (ljmin<=j .AND. j<=ljmax) THEN
             DO i = liminp1,limax
               IF ( h0_g(i,jh(i)) /= -valmanq) THEN
#if defined key_siggen || defined key_gencoord
                 DO k=1,kmax
                   dz3=f_dzu(h0(i,j),xe3(i,j),k,i,j)
                   phitemp(k,i,j)=phitemp(k,i,j)-halfdt*obc_coefrel*(temp(k,i,j)-obc_tn(k,i,j))*dz3*fwetp(i,j)
                   phisal(k,i,j) =    0.0_rlg -halfdt*obc_coefrel*(sal(k,i,j)-obc_sn(k,i,j))*dz3*fwetp(i,j)
#ifdef key_substance
#ifdef key_zerograd
                   DO iv=1,nv_adv
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                      phicon(iv,k,i,j) = phicon(iv,k,i,j) -halfdt*obc_coefrel*(cv_wat(iv,k,i,j)-obc_cv_wat_n(iv,k,i,j))*dz3*fwetp(i,j)
                    ENDIF
                   ENDDO
#else
                   phicon(:,k,i,j) = phicon(:,k,i,j) -halfdt*obc_coefrel*(cv_wat(:,k,i,j)-obc_cv_wat_n(:,k,i,j))*dz3*fwetp(i,j)
#endif
#endif
                 END DO
#else
                 phitemp(:,i,j)=phitemp(:,i,j)-halfdt*obc_coefrel*(temp(:,i,j)-obc_tn(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                 phisal(:,i,j) =    0.0_rlg -halfdt*obc_coefrel*(sal(:,i,j)-obc_sn(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#ifdef key_substance
                   DO iv=1,nv_adv
#ifdef key_zerograd
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                      phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_n(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                    ENDIF
#else
                      phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_n(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#endif
                   ENDDO
#endif

#endif
               END IF
             END DO
           ENDIF_MPI
         END DO
!$OMP END DO
       END IF

       IF (l_obc_west) THEN
!$OMP DO SCHEDULE(STATIC)
         DO j = ljmin,ljmax
           DO i = imin+2,imin+obc_width
             IF_MPI (limin<=i .AND. i<=limax) THEN
               IF ( h0_g(ig(j),j) /= -valmanq) THEN
#if defined key_siggen || defined key_gencoord
                 DO k=1,kmax
                   dz3=f_dzu(h0(i,j),xe3(i,j),k,i,j)
                   phitemp(k,i,j)=phitemp(k,i,j)-halfdt*obc_coefrel*(temp(k,i,j)-obc_tw(k,i,j))*dz3*fwetp(i,j)
                   phisal(k,i,j) =    0.0_rlg   -halfdt*obc_coefrel*(sal(k,i,j)-obc_sw(k,i,j))*dz3*fwetp(i,j)
#ifdef key_substance
#ifdef key_zerograd
                   DO iv=1,nv_adv
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                      phicon(iv,k,i,j) = phicon(iv,k,i,j) -halfdt*obc_coefrel*(cv_wat(iv,k,i,j)-obc_cv_wat_w(iv,k,i,j))*dz3*fwetp(i,j)
                    ENDIF
                   ENDDO
#else
                   phicon(:,k,i,j) = phicon(:,k,i,j) -halfdt*obc_coefrel*(cv_wat(:,k,i,j)-obc_cv_wat_w(:,k,i,j))*dz3*fwetp(i,j)
#endif
#endif

                 END DO
#else
                 phitemp(:,i,j)=phitemp(:,i,j)-halfdt*obc_coefrel*(temp(:,i,j)-obc_tw(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                 phisal(:,i,j) =    0.0_rlg   -halfdt*obc_coefrel*(sal(:,i,j)-obc_sw(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#ifdef key_substance
                   DO iv=1,nv_adv
#ifdef key_zerograd
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                      phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_w(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                    ENDIF
#else
                      phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_w(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#endif
                   ENDDO
#endif
#endif
               END IF
             ENDIF_MPI
           END DO
         END DO
!$OMP END DO
       END IF

       IF (l_obc_east) THEN
!$OMP DO SCHEDULE(STATIC)
         DO j = ljmin,ljmax
           DO i = imax-obc_width+1,imax
             IF_MPI (limin<=i .AND. i<=limax) THEN
               IF ( h0_g(id(j),j) /= -valmanq) THEN
#if defined key_siggen || defined key_gencoord
                 DO k=1,kmax
                   dz3=f_dzu(h0(i,j),xe3(i,j),k,i,j)
                   phitemp(k,i,j)=phitemp(k,i,j)-halfdt*obc_coefrel*(temp(k,i,j)-obc_te(k,i,j))*dz3*fwetp(i,j)
                   phisal(k,i,j) =    0.0_rlg   -halfdt*obc_coefrel*(sal(k,i,j)-obc_se(k,i,j))*dz3*fwetp(i,j)
#ifdef key_substance
#ifdef key_zerograd
                   DO iv=1,nv_adv
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                      phicon(iv,k,i,j) = phicon(iv,k,i,j) -halfdt*obc_coefrel*(cv_wat(iv,k,i,j)-obc_cv_wat_e(iv,k,i,j))*dz3*fwetp(i,j)
                    ENDIF
                   ENDDO
#else
                   phicon(:,k,i,j) = phicon(:,k,i,j) -halfdt*obc_coefrel*(cv_wat(:,k,i,j)-obc_cv_wat_e(:,k,i,j))*dz3*fwetp(i,j)
#endif
#endif
                 END DO
#else
                 phitemp(:,i,j)=phitemp(:,i,j)-halfdt*obc_coefrel*(temp(:,i,j)-obc_te(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                 phisal(:,i,j) =    0.0_rlg   -halfdt*obc_coefrel*(sal(:,i,j)-obc_se(:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#ifdef key_substance
                   DO iv=1,nv_adv
#ifdef key_zerograd
                    IF (obc_cv_name(iv) /= 'zerograd') THEN
                      phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_e(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
                    ENDIF
#else
                      phicon(iv,:,i,j) = phicon(iv,:,i,j) -halfdt*obc_coefrel*(cv_wat(iv,:,i,j)-obc_cv_wat_e(iv,:,i,j))*h0(i,j)*dsigu(:)*fwetp(i,j)
#endif
                   ENDDO
#endif
#endif
               END IF
             ENDIF_MPI
           END DO
         END DO
!$OMP END DO
       END IF

     END IF  ! ends relaxation


   ELSE

   ! -------------------------------
   ! constant radiative flux
   ! no infra-red and turbulent flux
   ! -------------------------------

   ! Initialize phitemp to 0 at each river discharge
   !------------------------------------------------

!$OMP SINGLE
     DO k = 1,nbfl
       i=iriv(k)
       j=jriv(k)
       IF_MPI (ljmin<=j .AND. j<=ljmax) THEN
         IF_MPI (limin<=i .AND. i<=limax) THEN
           phitemp(:,i,j) = 0.0_rlg
         ENDIF_MPI
       ENDIF_MPI
     ENDDO
!$OMP END SINGLE

   ! Initialize phitemp to the radiative flux at the surface
   !--------------------------------------------------------

!$OMP DO SCHEDULE(RUNTIME)
     DO j=ljminp1,ljmax
       DO i = liminp1,limax
         IF (h0(i,j)+xe3(i,j) > h0fond) phitemp(kmax,i,j)=halfdt*sflx_solarcst/(chp*rhoref)*fwetp(i,j)
       END DO
     END DO
!$OMP END DO

   ! apport de chaleur par les fleuves
   !-----------------------------------

!$OMP SINGLE
     tjour=jjulien+iheure/24.0_rlg+iminu/1440.0_rlg+isec/86400.0_rlg
     DO k = 1,nbfl
       IF (l_tempriv(k)) THEN
         tempfle=tmfle(k)
       ELSE
         tempfle=tmfle(k)+amtfle(k)*REAL(COS(2.0_rlg*pi/365.256363_rlg*(tjour-220.0_rlg)),rsh)*jhemisph
       END IF
       i=iriv(k)
       j=jriv(k)
       IF_MPI (ljmin<=j .AND. j<=ljmax) THEN
         IF_MPI (limin<=i .AND. i<=limax) THEN
           phitemp(:,i,j) = phitemp(:,i,j)+halfdt*phiriver(k,:)*tempfle/surf(i,j)
           phieau(1:kmax,i,j)=phieau(1:kmax,i,j)+phiriver(k,1:kmax)
         ENDIF_MPI
       ENDIF_MPI
     END DO
!$OMP END SINGLE

   END IF

   ! Atmospherical fluxes for substances
#ifdef key_substance
   IF (l_subflxatm) THEN
!$OMP DO SCHEDULE(RUNTIME)
     DO j=ljmin,ljmax
       DO iv=1,nv_adv
         IF(l_subflxatm_xyt) THEN
         ! Atmospherical fluxes for substances vary with time et space
            phicon_atm(iv,limin:limax,j)=halfdt*surf(limin:limax,j)*fwetp(limin:limax,j) * &
                                       ( sub_flx_atm_xyt(iv,limin:limax,j) +                     &
                                        0.001_rsh*rain(limin:limax,j)*cv_rain(iv) )
         ELSE
            phicon_atm(iv,limin:limax,j)=halfdt*surf(limin:limax,j)*fwetp(limin:limax,j) * &
                                         ( sub_flx_atm(iv) +                     &
                                           0.001_rsh*rain(limin:limax,j)*cv_rain(iv) )
         ENDIF

         DO i=limin,limax
           immer_top = 0.0_rsh
           d3=h0(i,j)+xe3(i,j)

           IF (d3 >= sflx_sub_atm_depth ) THEN
             k = kmax
             DO WHILE (immer_top < sflx_sub_atm_depth .AND. k/=0)
! attention on veut epaisseur de la couche
               ! immersion from the sea surface. The water thickness is xe - z (z being negative)
               immer_bot = - (f_zw(h0(i,j),xe3(i,j),k-1,i,j)-xe3(i,j))
               dz = f_dzu(h0(i,j),xe3(i,j),k,i,j)
               IF (immer_bot < sflx_sub_atm_depth) THEN
                 phicon(iv,k,i,j)=phicon(iv,k,i,j) + &
                                  phicon_atm(iv,i,j)*dz/sflx_sub_atm_depth
               ELSE
                 phicon(iv,k,i,j)=phicon(iv,k,i,j) + &
                                  phicon_atm(iv,i,j)*(sflx_sub_atm_depth-immer_top)/sflx_sub_atm_depth
               END IF
               immer_top = immer_bot
               k = k-1
             END DO
           ELSE IF(d3 > h0fond) THEN
             phicon(iv,:,i,j)=phicon(iv,:,i,j) + phicon_atm(iv,i,j)/kmax
           ELSE
           ! ATTENTION : actual assumption :  
           !       no atmospheric input if no water in the cell
           !       no rain either if d3<h0fond (see  before)
              phicon_atm(iv,i,j)=0.0_rsh
          ENDIF
         END DO
       END DO
     END DO
!$OMP END DO

   END IF
#endif

!$OMP END PARALLEL

END SUBROUTINE sflx_surf
