mixtln.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
    SUBROUTINE mixtl_n (O, OR, PLAT, PFSOL,PFNSOL,PSFTEAU,PSFU,PSFV,PSEATEMP)
!     #######################################################################
!
!!****  *MIXTLN (1D MODEL)*  
!!
!!    PURPOSE
!!    -------
!     Oceanic 1D model in TKE Closure scheme
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!      
!!    REFERENCE
!!    ---------
!     Gaspar et al, 1990 : A simple eddy kinetic energy model for simulations of 
!     the oceanic vertical mixing : Tests at station Papa and Long-term upper ocean 
!     study site, JGR, 95,C9, 16,179--16,193. 
!!      
!!    AUTHOR
!!    ------
!!     C. Lebeaupin  *Météo-France* 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original     02/2008
!!     01/2012 : H. Giordani, P. Peyrille
!!                       Add FLAGS: 
!!                       LREL_TS: relaxation on T, S
!!                       LREL_CUR: damping on current
!!                       FOR LREL_CUR: implicit and explicit codinf is made
!!                       for conveniency, DTREL=Ucur term, DSREL=VCURterm
!!                Corrections:
!!                   coriolis terms in current equation, 
!!                   richardson nb in diapycnal mixing, 
!!                   remove threshold value for mixing tendency
!!    07/2012, P. Le Moigne : CMO1D phasing
!!    09/2016, C. Lebeaupin Brossier: XSEATEND and initialization
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_ocean_n, ONLY : ocean_t
USE modd_ocean_rel_n, ONLY : ocean_rel_t
!
USE modd_csts
USE modd_ocean_csts
USE modd_ocean_grid
!
USE modd_surf_par,   ONLY : xundef
!
! Module containing relaxation fields
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
!
TYPE(ocean_t), INTENT(INOUT) :: O
TYPE(ocean_rel_t), INTENT(INOUT) :: OR
!
REAL, DIMENSION(:), INTENT(IN) :: PLAT
!
REAL, DIMENSION(:)  ,INTENT(IN)       :: PFSOL   ! solar flux (W/m2)
REAL, DIMENSION(:)  ,INTENT(IN)       :: PFNSOL  ! non solar flux (W/m2)
REAL, DIMENSION(:)  ,INTENT(IN)       :: PSFTEAU ! fresh water flux(kg/m2/s)
REAL, DIMENSION(:)  ,INTENT(IN)       :: PSFU    ! zonal stress (Pa)
REAL, DIMENSION(:)  ,INTENT(IN)       :: PSFV    ! meridian stress (Pa)
!
REAL, DIMENSION(:)  ,INTENT(OUT)    :: PSEATEMP! sea surface temperature (K)
!*      0.2    declarations of local variables
!
 COMPLEX, DIMENSION(NOCKMIN:NOCKMAX,2) :: ZUC           ! vecteur vent en ecriture complexe
 COMPLEX, DIMENSION(NOCKMIN:NOCKMAX) :: ZAI, ZBI, ZCI     !matrices pour resolution numérique
 COMPLEX, DIMENSION(NOCKMIN:NOCKMAX) :: ZAU, ZBU, ZCU, ZYU !matrices pour resolution numérique
 COMPLEX, DIMENSION(NOCKMIN:NOCKMAX) :: ZOMU,ZWU
!
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZSEAT, ZSEAS, ZSEAE, ZSEAV, ZSEAU
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZSEAT_REL, ZSEAS_REL, ZSEAV_REL, ZSEAU_REL
!
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZA, ZB, ZC, ZA2, ZB2, ZC2, ZYT, ZYS, ZYE        !matrices pour resolution numérique
!
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ADVT,ADVS !advection horiz. temperature and salinity
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ADVU,ADVV !advection horiz. of current
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ADVE      !advection of turbulent kinetic energy
!
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZTDTREL  ! Tendancy derived from  relxation (K/s)
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZSDTREL  ! ---- salinity   ---------- (%/s)
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZUDTREL  ! Tendancy term derived from relaxation on U current (m/s2)
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZVDTREL  !  ------------------------------------    V        (m/s2)
!
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZLE,ZKMEL,ZKMELM
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZKMES,ZKMED,ZKMEWM,ZKMEWS
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZPTH,ZPDY
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZTENDE,ZDIFFV
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZOMT, ZOMS, ZOME !vector for matrix inversion
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZWT, ZWS, ZWE
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZU,ZV,ZT,ZS,ZE
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZZDRHO
REAL, DIMENSION(NOCKMIN:NOCKMAX) :: ZDTFSOL
!
REAL :: ZDTFNSOL
!
REAL :: ZSFU, ZSFV, ZFNSOL, ZFSOL, ZSFTEAU, ZLAT
REAL :: ZSEAHMO, ZSEATEMP
!
REAL :: ZSEUIL,ZEMIN,ZEMAX,ZTEST
!
REAL :: ZF, ZEWS
REAL :: ZALG
REAL :: ZEE,ZPOT,ZXLME,ZXLPE,ZXROD,ZAUX,ZXDL
REAL :: ZDU,ZDV,ZRICH,ZDRHODZ
REAL :: ZT1, ZT2, ZT3, ZS1, ZS2
!
INTEGER :: IUP,IBOT
INTEGER :: J,JJ,JPT,JIN,IKHML
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!----------------------------------------------------------------------------
!       1.     Initializations
!__________________________________________________________________________
!Take account of an horizontale advection

!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MIXTL_N',0,zhook_handle)
iup=nockmin+1
ibot=nockmax
!
zseuil=2.23e-3
zemin=1.e-3
zemax=0.1
!
zt1=13.5
zt2=-0.19494
zt3=-0.49038e-2
zs1=32.6
zs2=0.77475
!
zalg = xg / xrhosw
!
zuc(:,:) = (0.,0.)
zai(:) = (0.,0.)
zbi(:) = (0.,0.)
zci(:) = (0.,0.)
zau(:) = (0.,0.)
zbu(:) = (0.,0.)
zcu(:) = (0.,0.)
zyu(:) = (0.,0.)
zomu(:) = (0.,0.)
zwu(:) = (0.,0.)
!
zseat(:) = 0.
zseas(:) = 0.
zseav(:) = 0.
zseae(:) = 0.
zseau(:) = 0.
!
zseat_rel(:) = 0.
zseas_rel(:) = 0.
zseav_rel(:) = 0.
zseau_rel(:) = 0.
!
za(:) = 0.
zb(:) = 0.
zc(:) = 0.
za2(:) = 0.
zb2(:) = 0.
zc2(:) = 0.
zyt(:) = 0.
zys(:) = 0.
zye(:) = 0.
!
advt(:) = 0.
advs(:) = 0.
advu(:) = 0.
advv(:) = 0.
adve(:) = 0.
!
ztdtrel(:) = 0.
zsdtrel(:) = 0.
zudtrel(:) = 0.
zvdtrel(:) = 0.
!
zle(:) = 0.
zkmel(:) = 0.
zkmelm(:) = 0.
zkmes(:) = 0.
zkmed(:) = 0.
zkmewm(:) = 0.
zkmews(:) = 0.
zpth(:) = 0.
zpdy(:) = 0.
ztende(:) = 0.
zdiffv(:) = 0.
zomt(:) = 0.
zoms(:) = 0.
zome(:) = 0.
zwt(:) = 0.
zws(:) = 0.
zwe(:) = 0.
zzdrho(:) = 0.
zdtfsol(:) = 0.
zt(:) = 0.
zs(:) = 0.    
zu(:) = 0.
zv(:) = 0.
ze(:) = 0.
!
!-------------------------------------------------------------------------------
!
!iterations on grid points
DO jpt=1,SIZE(pfsol)
!
  ztest=0.
  ikhml=1
!
  !simplified variables inside this loop
  zlat    = plat(jpt)
  zfsol   = pfsol(jpt)
  zfnsol  = pfnsol(jpt)
  zsfteau = psfteau(jpt)
  zsfu    = psfu(jpt)
  zsfv    = psfv(jpt)
  zews    = sqrt(zsfu**2+zsfv**2)/xrhosw
  zf      = 4.*xpi*sin(zlat*xpi/180.)/86400.

  zseat(:) = o%XSEAT(jpt,:)
  zseas(:) = o%XSEAS(jpt,:)
  zseau(:) = o%XSEAU(jpt,:)
  zseav(:) = o%XSEAV(jpt,:)    
  zseae(:) = o%XSEAE(jpt,:)
  !
  zseau_rel(:) = or%XSEAU_REL(jpt,:)
  zseav_rel(:) = or%XSEAV_REL(jpt,:)
  zseat_rel(:) = or%XSEAT_REL(jpt,:)
  zseas_rel(:) = or%XSEAS_REL(jpt,:)
  !
  zseahmo = 0.
  DO j=iup-1,ibot
    IF (j>=iup .AND. zseae(j)>=(zemin*sqrt(2.))) zseahmo = zseahmo-xdz1(j)
  ENDDO
  o%XSEAHMO(jpt) = zseahmo
  !
  !precalculation of DRHO
  zu(:) = 0.
  zv(:) = 0.
  zt(:) = 0.
  zs(:) = 0.
  ze(:) = 0.
  advt(:) = 0.
  advs(:) = 0.
  advu(:) = 0.
  advv(:) = 0.
  adve(:) = 0.
  zzdrho(:) = (zseat(:)-zt1)*(zt2+zt3*(zseat(:)-zt1)) + zs2*(zseas(:)-zs1)
  zudtrel(:) = 0.
  zvdtrel(:) = 0.
  ztdtrel(:) = 0.
  zsdtrel(:) = 0.
  zdtfsol(:) = 0.
  !
  zdtfnsol=0.
!
! Control print
!IF (LREL_CUR)  WRITE(*,*) "WARNING :: Damping on current will be done"
!IF (LREL_TS)  WRITE(*,*) "WARNING :: Relaxation on T, S ocean will be done "
!IF (LDIAPYCNAL)  WRITE(*,*) "WARNING :: diapycnal mixing has been activated"
!IF (LFLX_CORR) WRITE(*,*) "WARNING :: ocean fluxes correctin has been activated"
!
!------------------------------------------------------------------------------
!
!       2.     Oceanic vertical mixing 
!              -----------------------
!!
!!       2.a    Diapycnal mixing
!!              ----------------
!!
  DO j=ibot-1,iup,-1
    zkmes(j) =0.
    zkmewm(j)=0.
    zkmews(j)=0.
    zkmed(j) =0.
    IF ((ztest==0.).AND.(zseae(j)>=zseuil)) THEN
      ikhml=j
      ztest=1.
    ENDIF
  ENDDO

  !! Modif PP - HG : flag diapycnal
  IF (or%LDIAPYCNAL) THEN
    !
    DO j=ikhml,ibot-1
      zdrhodz=(zzdrho(j)-zzdrho(j+1))/xdz1(j)
      zdu=zseau(j+1)-zseau(j)
      zdv=zseav(j+1)-zseav(j)
      IF((zdu*zdu+zdv*zdv).LE.1.e-7) THEN 
        zrich = 0.8
      ELSE
        zrich = -zalg*zdrhodz/(zdu**2+zdv**2)/xk4(j)
      ENDIF  
!coefficient de mélange aux ondes internes
      zkmewm(j)=1.e-3
      zkmews(j)=1.e-4
!coefficient de mélange du au cisaillement
      IF(zrich>0.7) THEN
        zkmes(j) = 0.
      ELSEIF(zrich>=0.) THEN
        zkmes(j) = 5.e-3*(1.-(zrich/0.7)*(zrich/0.7))**3
      ELSE
        zkmes(j) = 5.e-3
      ENDIF
    ENDDO
    !plm si ldiapycnal=F zkmes non modofie et zrich ne sert a rien !
    !plm ELSE
    !plm  ZRICH=-ZALG*ZDRHODZ / (ZDU**2+ZDV**2) / XK4(J) 
    !  
  ENDIF
!
!       2.b    Mixing length and coefficient 
!              -----------------------------
!
  DO j=iup,ibot-1

    zee=zseae(j)**2/zalg
    zxrod=(zzdrho(j)+ zzdrho(j+1))*0.5   

    zxlme=0.
    zpot=0.
    jin=j
    DO jj=j+1,ibot
      zaux=zpot + xdz2(jj)*(zzdrho(jj)-zxrod)
      IF (zaux<zee) THEN
        jin=jj
        zpot=zaux
        zxlme=zxlme+xdz2(jj)
      ENDIF
    ENDDO
    IF (jin==j) THEN
      zxlme=2.*(xz2(j)-xzhoc(j+1))/(zzdrho(j+1)-zxrod)/zalg
      zxlme=sqrt(zxlme)*zseae(j)
    ELSE
      IF (jin/=ibot) THEN
        zxdl=(zee-zpot)/(zzdrho(jin+1)-zxrod)
        zxlme=zxlme+zxdl
      ENDIF
    ENDIF

    zxlpe=0.
    zpot=0.
    jin=j
    DO jj=j-1,nockmin,-1
      zaux=zpot + xdz2(jj+1)*(zxrod-zzdrho(jj+1))
      IF (zaux<zee) THEN
        jin=jj
        zpot=zaux
        zxlpe=zxlpe+xdz2(jj+1)
      ENDIF
    ENDDO
    IF (jin==j) THEN
      zxlpe=2.*(xz2(j)-xzhoc(j))/(zzdrho(j)-zxrod)/zalg
      zxlpe=sqrt(zxlpe)*zseae(j)
    ELSE 
      IF (jin/=nockmin) THEN
        zxdl=- (zee-zpot)/(zzdrho(jin)-zxrod)
        zxlpe=zxlpe+zxdl
      ENDIF
    ENDIF

    zle(j)=sqrt(zxlme*zxlpe)
    zkmel(j)=xckl*zle(j)*zseae(j)
  
  ENDDO

  zle(ibot)=zle(ibot-1)
  zkmel(ibot)=xckl*zle(ibot)*zseae(ibot)

  !first coef at all levels: because needed at j+1 further
  zkmelm(iup)=zkmel(iup)
  DO j=iup+1,ibot
    zkmelm(j)= (zkmel(j)+zkmel(j-1))/2.
  ENDDO

!--------------------------------------------------------------------------------------------------------
!
!!       2.c    Numerical resolution of evolution equations
!!              -------------------------------------------
!
  IF (or%LREL_CUR) THEN
    DO j=iup,ibot
      zudtrel(j) =  - (zseau(j)-zseau_rel(j))  / or%XTAU_REL 
      zvdtrel(j) =  - (zseav(j)-zseav_rel(j))  / or%XTAU_REL 
    ENDDO
  ENDIF
  !
  DO j=iup,ibot
    ! flux solaire
    zdtfsol(j) = xray(j)*zfsol/xdz2(j) 
  ENDDO
!
  IF (or%LREL_TS) THEN
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! RELAXATION IS MADE INSTEAD OF FLUX CORRECTION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
    DO j=iup,ibot
      ! flux non solaire
      ztdtrel(j) =  - (zseat(j)-zseat_rel(j)) / or%XTAU_REL
      zsdtrel(j) =  - (zseas(j)-zseas_rel(j)) / or%XTAU_REL
    ENDDO
  ENDIF
!
  IF (or%LFLX_CORR) THEN
    ! NO relaxation
    ! Barnier correction on surface fluxes          
    ! flux non solaire corrige
    ztdtrel(iup) =  or%XQCORR *(zseat_rel(iup)-zseat(iup)) / (xrhosw*xcpsw)
    zfnsol = zfnsol + ztdtrel(iup)
  ENDIF
!
! flux non solaire
  zdtfnsol = zfnsol/xdz2(iup) 
!
!!
!loop on levels
  za(:) = 0.
  zai(:) = 0.
  za2(:) = 0.
  zc(:) = 0.
  zci(:) = 0.
  zc2(:) = 0.
  DO j=iup,ibot

    IF (j<ibot) THEN
      IF (j>iup) THEN
        za(j) = o%XOCEAN_TSTEP * xk1(j) * (zkmel(j-1) + zkmes(j-1) + zkmews(j-1) + zkmed(j-1)) 
        zai(j) = o%XOCEAN_TSTEP * xk1(j) * (zkmel(j-1) + zkmes(j-1) + zkmewm(j-1) + zkmed(j-1))
        za2(j) = o%XOCEAN_TSTEP * xk2(j) * zkmelm(j)
      ENDIF
      zc(j) = o%XOCEAN_TSTEP * xk2(j) * (zkmel(j)   + zkmes(j)   + zkmews(j)   + zkmed(j)) 
      zci(j) = o%XOCEAN_TSTEP * xk2(j) * (zkmel(j)   + zkmes(j)   + zkmewm(j)   + zkmed(j))
      zc2(j) = o%XOCEAN_TSTEP * xk3(j) * zkmelm(j+1)
    ENDIF
  
    zb(j) = 1. - za(j) - zc(j)
    zbi(j) = 1. - zai(j) - zci(j)
    zb2(j) = 1. - za2(j) - zc2(j) + o%XOCEAN_TSTEP * zseae(j)/zle(j)/xzce

    zau(j)    = zai(j) * (1.,0.) 
    zcu(j)    = zci(j) * (1.,0.) 
    zbu(j)    = zbi(j) * (1.,0.) + o%XOCEAN_TSTEP * zf * xgama * (0.,1.)
    

    zomt(j)   = 1./(zb(j) - za(j) * zomt(j-1) * zc(j-1)) 
    zoms(j)   = 1./(zb(j) - za(j) * zoms(j-1) * zc(j-1))    
    zomu(j)   = 1./(zbu(j) - zau(j) * zomu(j-1) * zcu(j-1)) 
    zome(j)   = 1./(zb2(j) - za2(j) * zome(j-1) * zc2(j-1))
 

    zyt(j)    = zseat(j) + o%XOCEAN_TSTEP * (zdtfsol(j) + advt(j)) 
    zys(j)    = zseas(j) + o%XOCEAN_TSTEP * (             advs(j))  
    IF (or%LREL_TS) THEN
      IF (.NOT.or%LFLX_CORR) zyt(j) = zyt(j) + o%XOCEAN_TSTEP * ztdtrel(j)
      zys(j) = zys(j) + o%XOCEAN_TSTEP * zsdtrel(j)
    ENDIF
    
    zuc(j,1)  = zseau(j)*(1.,0.) + zseav(j)*(0.,1.)
    zyu(j)    = zuc(j,1) + o%XOCEAN_TSTEP * (zuc(j,1)*zf*(1.-xgama)*(0.,-1.) + advu(j)*(1.,0.) + advv(j)*(0.,1.))
    ! damping on current if LREL_CUR=T in explicit scheme
    IF (or%LREL_CUR) zyu(j) = zyu(j) + o%XOCEAN_TSTEP * (zudtrel(j)*(1.,0.) + zvdtrel(j)*(0.,1.))  

    IF (j<ibot) THEN
      zdrhodz   = (zzdrho(j)-zzdrho(j+1))/xdz1(j)
      zdu       = zseau(j+1)-zseau(j)
      zdv       = zseav(j+1)-zseav(j)

      zpth(j)   = xckl * zle(j) * zalg * zdrhodz
      zpdy(j)   = xckl * zle(j) * xk4(j) * (zdu**2+zdv**2)
    ELSE
      zpth(j) = zpth(j-1)
      zpdy(j) = zpdy(j-1)
    ENDIF

    zye(j)    = zseae(j) + o%XOCEAN_TSTEP * (0.5 * zseae(j)**2/zle(j)/xzce + adve(j)) + zpth(j) + zpdy(j)

    IF (j==iup) THEN
      zyt(j) = zyt(j) + o%XOCEAN_TSTEP * zdtfnsol
      zys(j) = zys(j) + o%XOCEAN_TSTEP * zseas(iup) * zsfteau / xdz2(iup)
      zyu(j) = zyu(j) - o%XOCEAN_TSTEP * ( zsfu*(1.,0.) + zsfv*(0.,1.) ) / xdz2(iup) / xrhosw
      zye(j) = zye(j) + o%XOCEAN_TSTEP * zews / xdz1(iup)
    ENDIF


    zwt(j)    = zomt(j) * (zyt(j)- za(j)*zwt(j-1))
    zws(j)    = zoms(j) * (zys(j)- za(j)*zws(j-1))
    zwu(j)    = zomu(j) * (zyu(j)- zau(j)*zwu(j-1))
    zwe(j)    = zome(j) * (zye(j)- za2(j)*zwe(j-1))

  ENDDO
!
!---------------------------------------------------------------------------------------------------
!
  zt(ibot)    = zwt(ibot)
  zs(ibot)    = zws(ibot)  
  zuc(ibot,2) = zwu(ibot)
  ze(ibot)    = zwe(ibot)  
  DO j=ibot-1,iup,-1
    zt(j)   = zwt(j) - zc(j) * zomt(j) * zt(j+1)
    zs(j)   = zws(j) - zc(j) * zoms(j) * zs(j+1)
    zuc(j,2) = zwu(j) - zcu(j) * zomu(j) * zuc(j+1,2)
    ze(j)   = zwe(j) - zc2(j) * zome(j) * ze(j+1)
  ENDDO
!
!
  DO j=iup,ibot
    zu(j)  = real(zuc(j,2))
    zv(j)  = aimag(zuc(j,2))
    !
    ze(j)  = max(zemin,ze(j))
  ! Transformation to preserve E <EMAX; secure if mixt crash
    ze(j)  = min(ze(j),zemax)
  !bilan TKE
    !ZTENDE(J) = (ZE(J)*ZE(J)-ZSEAE(J)**2)/O%XOCEAN_TSTEP
    !ZDIFFV(J) = ZTENDE(J) - ZSEAE(J)*(ZPDY(J) + ZPTH(J))
    !
    zseat(j)  = zt(j)
    zseas(j)  = zs(j)    
    zseau(j)  = zu(j)
    zseav(j)  = zv(j)
    zseae(j)  = ze(j)
  ENDDO
!
!------------------------------------------------------------------------------
!!       3.     New oceanic profiles
!!              --------------------
!!
  IF (o%LPROGSST) o%XSEATEND(jpt) = (zt(iup)-o%XSEAT(jpt,iup)) / o%XOCEAN_TSTEP
  zseat(nockmin)  = zt(iup)
  zseas(nockmin)  = zs(iup)  
  zseau(nockmin)  = zu(iup)
  zseav(nockmin)  = zv(iup)
  zseae(nockmin)  = ze(iup)

  !bathymetrie
  DO j=iup,ibot
    IF (o%XSEABATH(jpt,j)==0.) THEN
      zseat(j)  = zseat(j-1) 
      zseas(j)  = zseas(j-1)            
      zseau(j)  = zseau(j-1)
      zseav(j)  = zseav(j-1)
      zseae(j)  = zseae(j-1)
    ENDIF
  ENDDO
!
!SST diagnosticed with 1D oceanic model
  zseatemp = zseat(iup) + xtt
  pseatemp(jpt) = zseatemp
!
  DO j=iup,ibot
    o%XLE   (jpt,j) = zle(j)
    o%XLK   (jpt,j) = zle(j)
    o%XKMEL (jpt,j) = zkmel(j)
    o%XKMELM(jpt,j) = zkmelm(j)
  ENDDO

  DO j=iup-1,ibot
    o%XSEAT(jpt,j) = zseat(j)
    o%XSEAS(jpt,j) = zseas(j)
    o%XSEAU(jpt,j) = zseau(j)
    o%XSEAV(jpt,j) = zseav(j)
    o%XSEAE(jpt,j) = zseae(j)
    o%XDTFSOL(jpt,j) = zdtfsol(j)
  ENDDO

  o%XDTFNSOL(jpt) = zdtfnsol

ENDDO
!  
IF (lhook) CALL dr_hook('MIXTL_N',1,zhook_handle)
!
!!-------------------------------------------------------------------------------
!!
!ENDDO !end of iterations on sea surfex grid points
!!------------------------------------------------------------------------------
!
END SUBROUTINE mixtl_n