MODULE slab_heat_transp_mod

  ! Slab ocean : temperature tendencies due to horizontal diffusion
  ! and / or Ekman transport

  USE lmdz_grid_phy, ONLY: nbp_lon, nbp_lat, klon_glo
  USE lmdz_abort_physic, ONLY: abort_physic
  IMPLICIT NONE

  ! Variables copied over from dyn3d dynamics:
  REAL, SAVE, ALLOCATABLE :: fext(:) ! Coriolis f times cell area
  !$OMP THREADPRIVATE(fext)
  REAL, SAVE, ALLOCATABLE :: beta(:) ! df/dy
  !$OMP THREADPRIVATE(beta)
  REAL, SAVE, ALLOCATABLE :: unsairez(:) ! 1/cell area
  !$OMP THREADPRIVATE(unsairez)
  REAL, SAVE, ALLOCATABLE :: unsaire(:)
  !$OMP THREADPRIVATE(unsaire)
  REAL, SAVE, ALLOCATABLE :: cu(:) ! cell longitude dim (m)
  !$OMP THREADPRIVATE(cu)
  REAL, SAVE, ALLOCATABLE :: cv(:) ! cell latitude dim (m)
  !$OMP THREADPRIVATE(cv)
  REAL, SAVE, ALLOCATABLE :: cuvsurcv(:) ! cu/cv (v points)
  !$OMP THREADPRIVATE(cuvsurcv)
  REAL, SAVE, ALLOCATABLE :: cvusurcu(:) ! cv/cu (u points)
  !$OMP THREADPRIVATE(cvusurcu)
  REAL, SAVE, ALLOCATABLE :: aire(:) ! cell area
  !$OMP THREADPRIVATE(aire)
  REAL, SAVE :: apoln ! area of north pole points
  !$OMP THREADPRIVATE(apoln)
  REAL, SAVE :: apols ! area of south pole points
  !$OMP THREADPRIVATE(apols)
  REAL, SAVE, ALLOCATABLE :: aireu(:) ! area of u cells
  !$OMP THREADPRIVATE(aireu)
  REAL, SAVE, ALLOCATABLE :: airev(:) ! area of v cells
  !$OMP THREADPRIVATE(airev)

  ! Local parameters for slab transport
  LOGICAL, SAVE :: alpha_var ! variable coef for deep temp (1 layer)
  !$OMP THREADPRIVATE(alpha_var)
  LOGICAL, SAVE :: slab_upstream ! upstream scheme ? (1 layer)
  !$OMP THREADPRIVATE(slab_upstream)
  LOGICAL, SAVE :: slab_sverdrup ! use wind stress curl at equator
  !$OMP THREADPRIVATE(slab_sverdrup)
  LOGICAL, SAVE :: slab_tyeq ! use merid wind stress at equator
  !$OMP THREADPRIVATE(slab_tyeq)
  LOGICAL, SAVE :: ekman_zonadv ! use zonal advection by Ekman currents
  !$OMP THREADPRIVATE(ekman_zonadv)
  LOGICAL, SAVE :: ekman_zonavg ! zonally average wind stress
  !$OMP THREADPRIVATE(ekman_zonavg)

  REAL, SAVE :: alpham
  !$OMP THREADPRIVATE(alpham)
  REAL, SAVE :: gmkappa
  !$OMP THREADPRIVATE(gmkappa)
  REAL, SAVE :: gm_smax
  !$OMP THREADPRIVATE(gm_smax)

  ! geometry variables : f, beta, mask...
  REAL, SAVE, ALLOCATABLE :: zmasqu(:) ! continent mask for zonal mass flux
  !$OMP THREADPRIVATE(zmasqu)
  REAL, SAVE, ALLOCATABLE :: zmasqv(:) ! continent mask for merid mass flux
  !$OMP THREADPRIVATE(zmasqv)
  REAL, SAVE, ALLOCATABLE :: unsfv(:) ! 1/f, v points
  !$OMP THREADPRIVATE(unsfv)
  REAL, SAVE, ALLOCATABLE :: unsbv(:) ! 1/beta
  !$OMP THREADPRIVATE(unsbv)
  REAL, SAVE, ALLOCATABLE :: unsev(:) ! 1/epsilon (drag)
  !$OMP THREADPRIVATE(unsev)
  REAL, SAVE, ALLOCATABLE :: unsfu(:) ! 1/F, u points
  !$OMP THREADPRIVATE(unsfu)
  REAL, SAVE, ALLOCATABLE :: unseu(:)
  !$OMP THREADPRIVATE(unseu)

  ! Routines from dyn3d, valid on global dynamics grid only:
  PRIVATE :: gr_fi_dyn, gr_dyn_fi ! to go between 1D nd 2D horiz grid
  PRIVATE :: gr_scal_v, gr_v_scal, gr_scal_u ! change on 2D grid U,V, T points
  PRIVATE :: grad, diverg

CONTAINS

  SUBROUTINE ini_slab_transp_geom(ip1jm, ip1jmp1, unsairez_, fext_, unsaire_, &
          cu_, cuvsurcv_, cv_, cvusurcu_, &
          aire_, apoln_, apols_, &
          aireu_, airev_, rlatv, rad, omeg)
    ! number of points in lon, lat
    IMPLICIT NONE
    ! Routine copies some geometry variables from the dynamical core
    ! see global vars for meaning
    INTEGER, INTENT(IN) :: ip1jm
    INTEGER, INTENT(IN) :: ip1jmp1
    REAL, INTENT(IN) :: unsairez_(ip1jm)
    REAL, INTENT(IN) :: fext_(ip1jm)
    REAL, INTENT(IN) :: unsaire_(ip1jmp1)
    REAL, INTENT(IN) :: cu_(ip1jmp1)
    REAL, INTENT(IN) :: cuvsurcv_(ip1jm)
    REAL, INTENT(IN) :: cv_(ip1jm)
    REAL, INTENT(IN) :: cvusurcu_(ip1jmp1)
    REAL, INTENT(IN) :: aire_(ip1jmp1)
    REAL, INTENT(IN) :: apoln_
    REAL, INTENT(IN) :: apols_
    REAL, INTENT(IN) :: aireu_(ip1jmp1)
    REAL, INTENT(IN) :: airev_(ip1jm)
    REAL, INTENT(IN) :: rlatv(nbp_lat - 1)
    REAL, INTENT(IN) :: rad
    REAL, INTENT(IN) :: omeg

    CHARACTER (len = 20) :: modname = 'slab_heat_transp'
    CHARACTER (len = 80) :: abort_message

    ! Sanity check on dimensions
    IF ((ip1jm/=((nbp_lon + 1) * (nbp_lat - 1))).OR. &
            (ip1jmp1/=((nbp_lon + 1) * nbp_lat))) THEN
      abort_message = "ini_slab_transp_geom Error: wrong array sizes"
      CALL abort_physic(modname, abort_message, 1)
    endif
    ! Allocations could be done only on master process/thread...
    allocate(unsairez(ip1jm))
    unsairez(:) = unsairez_(:)
    allocate(fext(ip1jm))
    fext(:) = fext_(:)
    allocate(unsaire(ip1jmp1))
    unsaire(:) = unsaire_(:)
    allocate(cu(ip1jmp1))
    cu(:) = cu_(:)
    allocate(cuvsurcv(ip1jm))
    cuvsurcv(:) = cuvsurcv_(:)
    allocate(cv(ip1jm))
    cv(:) = cv_(:)
    allocate(cvusurcu(ip1jmp1))
    cvusurcu(:) = cvusurcu_(:)
    allocate(aire(ip1jmp1))
    aire(:) = aire_(:)
    apoln = apoln_
    apols = apols_
    allocate(aireu(ip1jmp1))
    aireu(:) = aireu_(:)
    allocate(airev(ip1jm))
    airev(:) = airev_(:)
    allocate(beta(nbp_lat - 1))
    beta(:) = 2 * omeg * cos(rlatv(:)) / rad

  END SUBROUTINE ini_slab_transp_geom

  SUBROUTINE ini_slab_transp(zmasq)

    !    USE lmdz_ioipsl_getin_p, ONLY: getin_p
    USE IOIPSL, ONLY: getin
    IMPLICIT NONE

    REAL zmasq(klon_glo) ! ocean / continent mask, 1=continent
    REAL zmasq_2d((nbp_lon + 1) * nbp_lat)
    REAL ff((nbp_lon + 1) * (nbp_lat - 1)) ! Coriolis parameter
    REAL eps ! epsilon friction timescale (s-1)
    INTEGER :: slab_ekman
    INTEGER i
    INTEGER :: iim, iip1, jjp1, ip1jm, ip1jmp1

    ! Some definition for grid size
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1)
    ip1jmp1 = (nbp_lon + 1) * nbp_lat
    iim = nbp_lon
    iip1 = nbp_lon + 1
    jjp1 = nbp_lat
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1)
    ip1jmp1 = (nbp_lon + 1) * nbp_lat

    ! Options for Heat transport
    ! Alpha variable?
    alpha_var = .FALSE.
    CALL getin('slab_alphav', alpha_var)
    PRINT *, 'alpha variable', alpha_var
    !  centered ou upstream scheme for meridional transport
    slab_upstream = .FALSE.
    CALL getin('slab_upstream', slab_upstream)
    PRINT *, 'upstream slab scheme', slab_upstream
    ! Sverdrup balance at equator ?
    slab_sverdrup = .TRUE.
    CALL getin('slab_sverdrup', slab_sverdrup)
    PRINT *, 'Sverdrup balance', slab_sverdrup
    ! Use tauy for meridional flux at equator ?
    slab_tyeq = .TRUE.
    CALL getin('slab_tyeq', slab_tyeq)
    PRINT *, 'Tauy forcing at equator', slab_tyeq
    ! Use tauy for meridional flux at equator ?
    ekman_zonadv = .TRUE.
    CALL getin('slab_ekman_zonadv', ekman_zonadv)
    PRINT *, 'Use Ekman flow in zonal direction', ekman_zonadv
    ! Use tauy for meridional flux at equator ?
    ekman_zonavg = .FALSE.
    CALL getin('slab_ekman_zonavg', ekman_zonavg)
    PRINT *, 'Use zonally-averaged wind stress ?', ekman_zonavg
    ! Value of alpha
    alpham = 2. / 3.
    CALL getin('slab_alpha', alpham)
    PRINT *, 'slab_alpha', alpham
    ! GM k coefficient (m2/s) for 2-layers
    gmkappa = 1000.
    CALL getin('slab_gmkappa', gmkappa)
    PRINT *, 'slab_gmkappa', gmkappa
    ! GM k coefficient (m2/s) for 2-layers
    gm_smax = 2e-3
    CALL getin('slab_gm_smax', gm_smax)
    PRINT *, 'slab_gm_smax', gm_smax
    ! -----------------------------------------------------------
    ! Define ocean / continent mask (no flux into continent cell)
    allocate(zmasqu(ip1jmp1))
    allocate(zmasqv(ip1jm))
    zmasqu = 1.
    zmasqv = 1.

    ! convert mask to 2D grid
    CALL gr_fi_dyn(1, iip1, jjp1, zmasq, zmasq_2d)
    ! put flux mask to 0 at boundaries of continent cells
    DO i = 1, ip1jmp1 - 1
      IF (zmasq_2d(i)>1e-5 .OR. zmasq_2d(i + 1)>1e-5) THEN
        zmasqu(i) = 0.
      ENDIF
    END DO
    DO i = iip1, ip1jmp1, iip1
      zmasqu(i) = zmasqu(i - iim)
    END DO
    DO i = 1, ip1jm
      IF (zmasq_2d(i)>1e-5 .OR. zmasq_2d(i + iip1)>1e-5) THEN
        zmasqv(i) = 0.
      END IF
    END DO

    ! -----------------------------------------------------------
    ! Coriolis and friction for Ekman transport
    slab_ekman = 2
    CALL getin("slab_ekman", slab_ekman)
    IF (slab_ekman>0) THEN
      allocate(unsfv(ip1jm))
      allocate(unsev(ip1jm))
      allocate(unsfu(ip1jmp1))
      allocate(unseu(ip1jmp1))
      allocate(unsbv(ip1jm))

      eps = 1e-5 ! Drag
      CALL getin('slab_eps', eps)
      PRINT *, 'epsilon=', eps
      ff = fext * unsairez ! Coriolis
      ! coefs to convert tau_x, tau_y to Ekman mass fluxes
      ! on 2D grid v points
      ! Compute correction factor [0 1] near the equator (f<<eps)
      IF (slab_sverdrup) THEN
        ! New formulation, sharper near equator, when eps gives Rossby Radius
        DO i = 1, ip1jm
          unsev(i) = exp(-ff(i) * ff(i) / eps**2)
        ENDDO
      ELSE
        DO i = 1, ip1jm
          unsev(i) = eps**2 / (ff(i) * ff(i) + eps**2)
        ENDDO
      END IF ! slab_sverdrup
      ! 1/beta
      DO i = 1, jjp1 - 1
        unsbv((i - 1) * iip1 + 1:i * iip1) = unsev((i - 1) * iip1 + 1:i * iip1) / beta(i)
      END DO
      ! 1/f
      ff = SIGN(MAX(ABS(ff), eps / 100.), ff) ! avoid value 0 at equator...
      DO i = 1, ip1jm
        unsfv(i) = (1. - unsev(i)) / ff(i)
      END DO
      ! compute values on 2D u grid
      ! 1/eps
      unsev(:) = unsev(:) / eps
      CALL gr_v_scal(1, unsfv, unsfu)
      CALL gr_v_scal(1, unsev, unseu)
    END IF

  END SUBROUTINE ini_slab_transp

  SUBROUTINE divgrad_phy(nlevs, temp, delta)
    ! Computes temperature tendency due to horizontal diffusion :
    ! T Laplacian, later multiplied by diffusion coef and time-step

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: nlevs ! nlevs : slab layers
    REAL, INTENT(IN) :: temp(klon_glo, nlevs) ! slab temperature
    REAL, INTENT(OUT) :: delta(klon_glo, nlevs) ! temp laplacian (heat flux div.)
    REAL :: delta_2d((nbp_lon + 1) * nbp_lat, nlevs)
    REAL ghx((nbp_lon + 1) * nbp_lat, nlevs), ghy((nbp_lon + 1) * (nbp_lat - 1), nlevs)
    INTEGER :: ll, iip1, jjp1

    iip1 = nbp_lon + 1
    jjp1 = nbp_lat

    ! transpose temp to 2D horiz. grid
    CALL gr_fi_dyn(nlevs, iip1, jjp1, temp, delta_2d)
    ! computes gradient (proportional to heat flx)
    CALL grad(nlevs, delta_2d, ghx, ghy)
    ! put flux to 0 at ocean / continent boundary
    DO ll = 1, nlevs
      ghx(:, ll) = ghx(:, ll) * zmasqu
      ghy(:, ll) = ghy(:, ll) * zmasqv
    END DO
    ! flux divergence
    CALL diverg(nlevs, ghx, ghy, delta_2d)
    ! laplacian back to 1D grid
    CALL gr_dyn_fi(nlevs, iip1, jjp1, delta_2d, delta)

  END SUBROUTINE divgrad_phy

  SUBROUTINE slab_ekman1(tx_phy, ty_phy, ts_phy, dt_phy)
    ! 1.5 Layer Ekman transport temperature tendency
    ! note : tendency dt later multiplied by (delta t)/(rho.H)
    ! to convert from divergence of heat fluxes to T

    IMPLICIT NONE

    ! tx, ty : wind stress (different grids)
    ! fluxm, fluz : mass *or heat* fluxes
    ! dt : temperature tendency
    INTEGER ij

    ! ts surface temp, td deep temp (diagnosed)
    REAL ts_phy(klon_glo)
    REAL ts((nbp_lon + 1) * nbp_lat), td((nbp_lon + 1) * nbp_lat)
    ! zonal and meridional wind stress
    REAL tx_phy(klon_glo), ty_phy(klon_glo)
    REAL tyu((nbp_lon + 1) * nbp_lat), txu((nbp_lon + 1) * nbp_lat)
    REAL txv((nbp_lon + 1) * (nbp_lat - 1)), tyv((nbp_lon + 1) * (nbp_lat - 1))
    REAL tcurl((nbp_lon + 1) * (nbp_lat - 1))
    ! zonal and meridional Ekman mass fluxes at u, v points (2D grid)
    REAL fluxz((nbp_lon + 1) * nbp_lat), fluxm((nbp_lon + 1) * (nbp_lat - 1))
    ! vertical  and absolute mass fluxes (to estimate alpha)
    REAL fluxv((nbp_lon + 1) * nbp_lat), fluxt((nbp_lon + 1) * (nbp_lat - 1))
    ! temperature tendency
    REAL dt((nbp_lon + 1) * nbp_lat), dt_phy(klon_glo)
    REAL alpha((nbp_lon + 1) * nbp_lat) ! deep temperature coef

    INTEGER iim, iip1, iip2, jjp1, ip1jm, ip1jmi1, ip1jmp1

    ! Grid definitions
    iim = nbp_lon
    iip1 = nbp_lon + 1
    iip2 = nbp_lon + 2
    jjp1 = nbp_lat
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1

    ! Convert taux,y to 2D  scalar grid
    ! Note: 2D grid size = iim*jjm. iip1=iim+1
    ! First and last points in zonal direction are the same
    ! we use 1 index ij from 1 to (iim+1)*(jjm+1)
    ! north and south poles
    tx_phy(1) = 0.
    tx_phy(klon_glo) = 0.
    ty_phy(1) = 0.
    ty_phy(klon_glo) = 0.
    CALL gr_fi_dyn(1, iip1, jjp1, tx_phy, txu)
    CALL gr_fi_dyn(1, iip1, jjp1, ty_phy, tyu)
    ! convert to u,v grid (Arakawa C)
    ! Multiply by f or eps to get mass flux
    ! Meridional mass flux
    CALL gr_scal_v(1, txu, txv) ! wind stress at v points
    IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
      tcurl = (txu(1:ip1jm) - txu(iip2:ip1jmp1)) / cv(:)
      fluxm = -tcurl * unsbv - txv * unsfv ! in kg.s-1.m-1 (zonal distance)
    ELSE
      CALL gr_scal_v(1, tyu, tyv)
      fluxm = tyv * unsev - txv * unsfv ! in kg.s-1.m-1 (zonal distance)
    ENDIF
    ! Zonal mass flux
    CALL gr_scal_u(1, txu, txu) ! wind stress at u points
    CALL gr_scal_u(1, tyu, tyu)
    fluxz = tyu * unsfu + txu * unseu

    ! Correct flux: continent mask and horiz grid size
    ! multiply m-flux by mask and dx: flux in kg.s-1
    fluxm = fluxm * cv * cuvsurcv * zmasqv
    ! multiply z-flux by mask and dy: flux in kg.s-1
    fluxz = fluxz * cu * cvusurcu * zmasqu

    ! Compute vertical  and absolute mass flux (for variable alpha)
    IF (alpha_var) THEN
      DO ij = iip2, ip1jm
        fluxv(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
        fluxt(ij) = ABS(fluxz(ij)) + ABS(fluxz(ij - 1)) &
                + ABS(fluxm(ij)) + ABS(fluxm(ij - iip1))
      ENDDO
      DO ij = iip1, ip1jmi1, iip1
        fluxt(ij + 1) = fluxt(ij + iip1)
        fluxv(ij + 1) = fluxv(ij + iip1)
      END DO
      fluxt(1) = SUM(ABS(fluxm(1:iim)))
      fluxt(ip1jmp1) = SUM(ABS(fluxm(ip1jm - iim:ip1jm - 1)))
      fluxv(1) = -SUM(fluxm(1:iim))
      fluxv(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))
      fluxt = MAX(fluxt, 1.e-10)
    ENDIF

    ! Compute alpha coefficient.
    ! Tdeep = Tsurf * alpha + 271.35 * (1-alpha)
    IF (alpha_var) THEN
      ! increase alpha (and Tdeep) in downwelling regions
      ! and decrease in upwelling regions
      ! to avoid "hot spots" where there is surface convergence
      DO ij = iip2, ip1jm
        alpha(ij) = alpham - fluxv(ij) / fluxt(ij) * (1. - alpham)
      ENDDO
      alpha(1:iip1) = alpham - fluxv(1) / fluxt(1) * (1. - alpham)
      alpha(ip1jm + 1:ip1jmp1) = alpham - fluxv(ip1jmp1) / fluxt(ip1jmp1) * (1. - alpham)
    ELSE
      alpha(:) = alpham
      ! Tsurf-Tdeep ~ 10° in the Tropics
    ENDIF

    ! Estimate deep temperature
    CALL gr_fi_dyn(1, iip1, jjp1, ts_phy, ts)
    DO ij = 1, ip1jmp1
      td(ij) = 271.35 + (ts(ij) - 271.35) * alpha(ij)
      td(ij) = MIN(td(ij), ts(ij))
    END DO

    ! Meridional heat flux: multiply mass flux by (ts-td)
    ! flux in K.kg.s-1
    IF (slab_upstream) THEN
      ! upstream scheme to avoid hot spots
      DO ij = 1, ip1jm
        IF (fluxm(ij)>=0.) THEN
          fluxm(ij) = fluxm(ij) * (ts(ij + iip1) - td(ij))
        ELSE
          fluxm(ij) = fluxm(ij) * (ts(ij) - td(ij + iip1))
        END IF
      END DO
    ELSE
      ! centered scheme better in mid-latitudes
      DO ij = 1, ip1jm
        fluxm(ij) = fluxm(ij) * (ts(ij + iip1) + ts(ij) - td(ij) - td(ij + iip1)) / 2.
      END DO
    ENDIF

    ! Zonal heat flux
    ! upstream scheme
    DO ij = iip2, ip1jm
      fluxz(ij) = fluxz(ij) * (ts(ij) + ts(ij + 1) - td(ij + 1) - td(ij)) / 2.
    END DO
    DO ij = iip1 * 2, ip1jmp1, iip1
      fluxz(ij) = fluxz(ij - iim)
    END DO

    ! temperature tendency = divergence of heat fluxes
    ! dt in K.s-1.kg.m-2 (T trend times mass/horiz surface)
    DO ij = iip2, ip1jm
      dt(ij) = (fluxz(ij - 1) - fluxz(ij) + fluxm(ij) - fluxm(ij - iip1)) &
              / aire(ij) ! aire : grid area
    END DO
    DO ij = iip1, ip1jmi1, iip1
      dt(ij + 1) = dt(ij + iip1)
    END DO
    ! special treatment at the Poles
    dt(1) = SUM(fluxm(1:iim)) / apoln
    dt(ip1jmp1) = -SUM(fluxm(ip1jm - iim:ip1jm - 1)) / apols
    dt(2:iip1) = dt(1)
    dt(ip1jm + 1:ip1jmp1) = dt(ip1jmp1)

    ! tendencies back to 1D grid
    CALL gr_dyn_fi(1, iip1, jjp1, dt, dt_phy)

  END SUBROUTINE slab_ekman1

  SUBROUTINE slab_ekman2(tx_phy, ty_phy, ts_phy, dt_phy_ek, dt_phy_gm, slab_gm)
    ! Temperature tendency for 2-layers slab ocean
    ! note : tendency dt later multiplied by (delta time)/(rho.H)
    ! to convert from divergence of heat fluxes to T

    ! 11/16 : Inclusion of GM-like eddy advection

    IMPLICIT NONE

    LOGICAL, INTENT(IN) :: slab_gm
    ! Here, temperature and flux variables are on 2 layers
    INTEGER ij

    ! wind stress variables
    REAL tx_phy(klon_glo), ty_phy(klon_glo)
    REAL txv((nbp_lon + 1) * (nbp_lat - 1)), tyv((nbp_lon + 1) * (nbp_lat - 1))
    REAL tyu((nbp_lon + 1) * nbp_lat), txu((nbp_lon + 1) * nbp_lat)
    REAL tcurl((nbp_lon + 1) * (nbp_lat - 1))
    ! slab temperature on  1D, 2D grid
    REAL ts_phy(klon_glo, 2), ts((nbp_lon + 1) * nbp_lat, 2)
    ! Temperature gradient, v-points
    REAL dty((nbp_lon + 1) * (nbp_lat - 1)), dtx((nbp_lon + 1) * nbp_lat)
    ! Vertical temperature difference, V-points
    REAL dtz((nbp_lon + 1) * (nbp_lat - 1))
    ! zonal and meridional mass fluxes at u, v points (2D grid)
    REAL fluxz((nbp_lon + 1) * nbp_lat), fluxm((nbp_lon + 1) * (nbp_lat - 1))
    ! vertical mass flux between the 2 layers
    REAL fluxv_ek((nbp_lon + 1) * nbp_lat)
    REAL fluxv_gm((nbp_lon + 1) * nbp_lat)
    ! zonal and meridional heat fluxes
    REAL fluxtz((nbp_lon + 1) * nbp_lat, 2)
    REAL fluxtm((nbp_lon + 1) * (nbp_lat - 1), 2)
    ! temperature tendency (in K.s-1.kg.m-2)
    REAL dt_ek((nbp_lon + 1) * nbp_lat, 2), dt_phy_ek(klon_glo, 2)
    REAL dt_gm((nbp_lon + 1) * nbp_lat, 2), dt_phy_gm(klon_glo, 2)
    ! helper vars
    REAL zonavg, fluxv
    REAL, PARAMETER :: sea_den = 1025. ! sea water density

    INTEGER iim, iip1, iip2, jjp1, ip1jm, ip1jmi1, ip1jmp1

    ! Grid definitions
    iim = nbp_lon
    iip1 = nbp_lon + 1
    iip2 = nbp_lon + 2
    jjp1 = nbp_lat
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
    ! Convert temperature to 2D grid
    CALL gr_fi_dyn(2, iip1, jjp1, ts_phy, ts)

    ! ------------------------------------
    ! Ekman mass fluxes and Temp tendency
    ! ------------------------------------
    ! Convert taux,y to 2D  scalar grid
    ! north and south poles tx,ty no meaning
    tx_phy(1) = 0.
    tx_phy(klon_glo) = 0.
    ty_phy(1) = 0.
    ty_phy(klon_glo) = 0.
    CALL gr_fi_dyn(1, iip1, jjp1, tx_phy, txu)
    CALL gr_fi_dyn(1, iip1, jjp1, ty_phy, tyu)
    IF (ekman_zonavg) THEN ! use zonal average of wind stress
      DO ij = 1, jjp1 - 2
        zonavg = SUM(txu(ij * iip1 + 1:ij * iip1 + iim)) / iim
        txu(ij * iip1 + 1:(ij + 1) * iip1) = zonavg
        zonavg = SUM(tyu(ij * iip1 + 1:ij * iip1 + iim)) / iim
        tyu(ij * iip1 + 1:(ij + 1) * iip1) = zonavg
      END DO
    END IF

    ! Divide taux,y by f or eps, and convert to 2D u,v grids
    ! (Arakawa C grid)
    ! Meridional flux
    CALL gr_scal_v(1, txu, txv) ! wind stress at v points
    fluxm = -txv * unsfv ! in kg.s-1.m-1 (zonal distance)
    IF (slab_sverdrup) THEN ! Sverdrup bal. near equator
      tcurl = (txu(1:ip1jm) - txu(iip2:ip1jmp1)) / cv(:) ! dtx/dy
      !poles curl = 0
      tcurl(1:iip1) = 0.
      tcurl(ip1jmi1 + 1:ip1jm) = 0.
      fluxm = fluxm - tcurl * unsbv
    ENDIF
    IF (slab_tyeq) THEN ! meridional wind forcing at equator
      CALL gr_scal_v(1, tyu, tyv)
      fluxm = fluxm + tyv * unsev ! in kg.s-1.m-1 (zonal distance)
    ENDIF
    !  apply continent mask, multiply by horiz grid dimension
    fluxm = fluxm * cv * cuvsurcv * zmasqv

    ! Zonal flux
    IF (ekman_zonadv) THEN
      CALL gr_scal_u(1, txu, txu) ! wind stress at u points
      CALL gr_scal_u(1, tyu, tyu)
      fluxz = tyu * unsfu + txu * unseu
      !  apply continent mask, multiply by horiz grid dimension
      fluxz = fluxz * cu * cvusurcu * zmasqu
    END IF

    !  Vertical mass flux from mass budget (divergence of horiz fluxes)
    IF (ekman_zonadv) THEN
      DO ij = iip2, ip1jm
        fluxv_ek(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
      ENDDO
    ELSE
      DO ij = iip2, ip1jm
        fluxv_ek(ij) = -fluxm(ij) + fluxm(ij - iip1)
      ENDDO
    END IF
    DO ij = iip1, ip1jmi1, iip1
      fluxv_ek(ij + 1) = fluxv_ek(ij + iip1)
    END DO
    !  vertical mass flux at Poles
    fluxv_ek(1) = -SUM(fluxm(1:iim))
    fluxv_ek(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))

    ! Meridional heat fluxes
    DO ij = 1, ip1jm
      ! centered scheme
      fluxtm(ij, 1) = fluxm(ij) * (ts(ij + iip1, 1) + ts(ij, 1)) / 2.
      fluxtm(ij, 2) = -fluxm(ij) * (ts(ij + iip1, 2) + ts(ij, 2)) / 2.
    END DO

    ! Zonal heat fluxes
    ! Schema upstream
    IF (ekman_zonadv) THEN
      DO ij = iip2, ip1jm
        IF (fluxz(ij)>=0.) THEN
          fluxtz(ij, 1) = fluxz(ij) * ts(ij, 1)
          fluxtz(ij, 2) = -fluxz(ij) * ts(ij + 1, 2)
        ELSE
          fluxtz(ij, 1) = fluxz(ij) * ts(ij + 1, 1)
          fluxtz(ij, 2) = -fluxz(ij) * ts(ij, 2)
        ENDIF
      ENDDO
      DO ij = iip1 * 2, ip1jmp1, iip1
        fluxtz(ij, :) = fluxtz(ij - iim, :)
      END DO
    ELSE
      fluxtz(:, :) = 0.
    ENDIF

    ! Temperature tendency, horizontal advection:
    DO ij = iip2, ip1jm
      dt_ek(ij, :) = fluxtz(ij - 1, :) - fluxtz(ij, :) &
              + fluxtm(ij, :) - fluxtm(ij - iip1, :)
    END DO
    ! Poles
    dt_ek(1, :) = SUM(fluxtm(1:iim, :), dim = 1)
    dt_ek(ip1jmp1, :) = -SUM(fluxtm(ip1jm - iim:ip1jm - 1, :), dim = 1)

    ! ------------------------------------
    ! GM mass fluxes and Temp tendency
    ! ------------------------------------
    IF (slab_gm) THEN
      ! Vertical Temperature difference T1-T2 on v-grid points
      CALL gr_scal_v(1, ts(:, 1) - ts(:, 2), dtz)
      dtz(:) = MAX(dtz(:), 0.25)
      ! Horizontal Temperature differences
      CALL grad(1, (ts(:, 1) + ts(:, 2)) / 2., dtx, dty)
      ! Meridional flux = -k.s (s=dyT/dzT)
      ! Continent mask, multiply by dz/dy
      fluxm = dty / dtz * 500. * cuvsurcv * zmasqv
      ! slope limitation, multiply by kappa
      fluxm = -gmkappa * SIGN(MIN(ABS(fluxm), gm_smax * cv * cuvsurcv), dty)
      ! convert to kg/s
      fluxm(:) = fluxm(:) * sea_den
      ! Zonal flux = 0. (temporary)
      fluxz(:) = 0.
      !  Vertical mass flux from mass budget (divergence of horiz fluxes)
      DO ij = iip2, ip1jm
        fluxv_gm(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
      ENDDO
      DO ij = iip1, ip1jmi1, iip1
        fluxv_gm(ij + 1) = fluxv_gm(ij + iip1)
      END DO
      !  vertical mass flux at Poles
      fluxv_gm(1) = -SUM(fluxm(1:iim))
      fluxv_gm(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))

      ! Meridional heat fluxes
      DO ij = 1, ip1jm
        ! centered scheme
        fluxtm(ij, 1) = fluxm(ij) * (ts(ij + iip1, 1) + ts(ij, 1)) / 2.
        fluxtm(ij, 2) = -fluxm(ij) * (ts(ij + iip1, 2) + ts(ij, 2)) / 2.
      END DO

      ! Zonal heat fluxes
      ! Schema upstream
      DO ij = iip2, ip1jm
        IF (fluxz(ij)>=0.) THEN
          fluxtz(ij, 1) = fluxz(ij) * ts(ij, 1)
          fluxtz(ij, 2) = -fluxz(ij) * ts(ij + 1, 2)
        ELSE
          fluxtz(ij, 1) = fluxz(ij) * ts(ij + 1, 1)
          fluxtz(ij, 2) = -fluxz(ij) * ts(ij, 2)
        ENDIF
      ENDDO
      DO ij = iip1 * 2, ip1jmp1, iip1
        fluxtz(ij, :) = fluxtz(ij - iim, :)
      END DO

      ! Temperature tendency :
      ! divergence of horizontal heat fluxes
      DO ij = iip2, ip1jm
        dt_gm(ij, :) = fluxtz(ij - 1, :) - fluxtz(ij, :) &
                + fluxtm(ij, :) - fluxtm(ij - iip1, :)
      END DO
      ! Poles
      dt_gm(1, :) = SUM(fluxtm(1:iim, :), dim = 1)
      dt_gm(ip1jmp1, :) = -SUM(fluxtm(ip1jm - iim:ip1jm - 1, :), dim = 1)
    ELSE
      dt_gm(:, :) = 0.
      fluxv_gm(:) = 0.
    ENDIF ! slab_gm

    ! ------------------------------------
    ! Temp tendency from vertical advection
    ! Divide by cell area
    ! ------------------------------------
    ! vertical heat flux = mass flux * T, upstream scheme
    DO ij = iip2, ip1jm
      fluxv = fluxv_ek(ij) + fluxv_gm(ij) ! net flux, needed for upstream scheme
      IF (fluxv>0.) THEN
        dt_ek(ij, 1) = dt_ek(ij, 1) + fluxv_ek(ij) * ts(ij, 2)
        dt_ek(ij, 2) = dt_ek(ij, 2) - fluxv_ek(ij) * ts(ij, 2)
        dt_gm(ij, 1) = dt_gm(ij, 1) + fluxv_gm(ij) * ts(ij, 2)
        dt_gm(ij, 2) = dt_gm(ij, 2) - fluxv_gm(ij) * ts(ij, 2)
      ELSE
        dt_ek(ij, 1) = dt_ek(ij, 1) + fluxv_ek(ij) * ts(ij, 1)
        dt_ek(ij, 2) = dt_ek(ij, 2) - fluxv_ek(ij) * ts(ij, 1)
        dt_gm(ij, 1) = dt_gm(ij, 1) + fluxv_gm(ij) * ts(ij, 1)
        dt_gm(ij, 2) = dt_gm(ij, 2) - fluxv_gm(ij) * ts(ij, 1)
      ENDIF
      ! divide by cell area
      dt_ek(ij, :) = dt_ek(ij, :) / aire(ij)
      dt_gm(ij, :) = dt_gm(ij, :) / aire(ij)
    END DO
    ! North Pole
    fluxv = fluxv_ek(1) + fluxv_gm(1)
    IF (fluxv>0.) THEN
      dt_ek(1, 1) = dt_ek(1, 1) + fluxv_ek(1) * ts(1, 2)
      dt_ek(1, 2) = dt_ek(1, 2) - fluxv_ek(1) * ts(1, 2)
      dt_gm(1, 1) = dt_gm(1, 1) + fluxv_gm(1) * ts(1, 2)
      dt_gm(1, 2) = dt_gm(1, 2) - fluxv_gm(1) * ts(1, 2)
    ELSE
      dt_ek(1, 1) = dt_ek(1, 1) + fluxv_ek(1) * ts(1, 1)
      dt_ek(1, 2) = dt_ek(1, 2) - fluxv_ek(1) * ts(1, 1)
      dt_gm(1, 1) = dt_gm(1, 1) + fluxv_gm(1) * ts(1, 1)
      dt_gm(1, 2) = dt_gm(1, 2) - fluxv_gm(1) * ts(1, 1)
    ENDIF
    dt_ek(1, :) = dt_ek(1, :) / apoln
    dt_gm(1, :) = dt_gm(1, :) / apoln
    ! South pole
    fluxv = fluxv_ek(ip1jmp1) + fluxv_gm(ip1jmp1)
    IF (fluxv>0.) THEN
      dt_ek(ip1jmp1, 1) = dt_ek(ip1jmp1, 1) + fluxv_ek(ip1jmp1) * ts(ip1jmp1, 2)
      dt_ek(ip1jmp1, 2) = dt_ek(ip1jmp1, 2) - fluxv_ek(ip1jmp1) * ts(ip1jmp1, 2)
      dt_gm(ip1jmp1, 1) = dt_gm(ip1jmp1, 1) + fluxv_gm(ip1jmp1) * ts(ip1jmp1, 2)
      dt_gm(ip1jmp1, 2) = dt_gm(ip1jmp1, 2) - fluxv_gm(ip1jmp1) * ts(ip1jmp1, 2)
    ELSE
      dt_ek(ip1jmp1, 1) = dt_ek(ip1jmp1, 1) + fluxv_ek(ip1jmp1) * ts(ip1jmp1, 1)
      dt_ek(ip1jmp1, 2) = dt_ek(ip1jmp1, 2) - fluxv_ek(ip1jmp1) * ts(ip1jmp1, 1)
      dt_gm(ip1jmp1, 1) = dt_gm(ip1jmp1, 1) + fluxv_gm(ip1jmp1) * ts(ip1jmp1, 1)
      dt_gm(ip1jmp1, 2) = dt_gm(ip1jmp1, 2) - fluxv_gm(ip1jmp1) * ts(ip1jmp1, 1)
    ENDIF
    dt_ek(ip1jmp1, :) = dt_ek(ip1jmp1, :) / apols
    dt_gm(ip1jmp1, :) = dt_gm(ip1jmp1, :) / apols

    dt_ek(2:iip1, 1) = dt_ek(1, 1)
    dt_ek(2:iip1, 2) = dt_ek(1, 2)
    dt_gm(2:iip1, 1) = dt_gm(1, 1)
    dt_gm(2:iip1, 2) = dt_gm(1, 2)
    dt_ek(ip1jm + 1:ip1jmp1, 1) = dt_ek(ip1jmp1, 1)
    dt_ek(ip1jm + 1:ip1jmp1, 2) = dt_ek(ip1jmp1, 2)
    dt_gm(ip1jm + 1:ip1jmp1, 1) = dt_gm(ip1jmp1, 1)
    dt_gm(ip1jm + 1:ip1jmp1, 2) = dt_gm(ip1jmp1, 2)

    DO ij = iip1, ip1jmi1, iip1
      dt_gm(ij + 1, :) = dt_gm(ij + iip1, :)
      dt_ek(ij + 1, :) = dt_ek(ij + iip1, :)
    END DO

    ! T tendency back to 1D grid...
    CALL gr_dyn_fi(2, iip1, jjp1, dt_ek, dt_phy_ek)
    CALL gr_dyn_fi(2, iip1, jjp1, dt_gm, dt_phy_gm)

  END SUBROUTINE slab_ekman2

  SUBROUTINE slab_gmdiff(ts_phy, dt_phy)
    ! Temperature tendency for 2-layers slab ocean
    ! Due to Gent-McWilliams type eddy-induced advection

    IMPLICIT NONE

    ! Here, temperature and flux variables are on 2 layers
    INTEGER ij
    ! Temperature gradient, v-points
    REAL dty((nbp_lon + 1) * (nbp_lat - 1)), dtx((nbp_lon + 1) * nbp_lat)
    ! Vertical temperature difference, V-points
    REAL dtz((nbp_lon + 1) * (nbp_lat - 1))
    ! slab temperature on  1D, 2D grid
    REAL ts_phy(klon_glo, 2), ts((nbp_lon + 1) * nbp_lat, 2)
    ! zonal and meridional mass fluxes at u, v points (2D grid)
    REAL fluxz((nbp_lon + 1) * nbp_lat), fluxm((nbp_lon + 1) * (nbp_lat - 1))
    ! vertical mass flux between the 2 layers
    REAL fluxv((nbp_lon + 1) * nbp_lat)
    ! zonal and meridional heat fluxes
    REAL fluxtz((nbp_lon + 1) * nbp_lat, 2)
    REAL fluxtm((nbp_lon + 1) * (nbp_lat - 1), 2)
    ! temperature tendency (in K.s-1.kg.m-2)
    REAL dt((nbp_lon + 1) * nbp_lat, 2), dt_phy(klon_glo, 2)

    INTEGER iim, iip1, iip2, jjp1, ip1jm, ip1jmi1, ip1jmp1

    ! Grid definitions
    iim = nbp_lon
    iip1 = nbp_lon + 1
    iip2 = nbp_lon + 2
    jjp1 = nbp_lat
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1

    ! Convert temperature to 2D grid
    CALL gr_fi_dyn(2, iip1, jjp1, ts_phy, ts)
    ! Vertical Temperature difference T1-T2 on v-grid points
    CALL gr_scal_v(1, ts(:, 1) - ts(:, 2), dtz)
    dtz(:) = MAX(dtz(:), 0.25)
    ! Horizontal Temperature differences
    CALL grad(1, (ts(:, 1) + ts(:, 2)) / 2., dtx, dty)
    ! Meridional flux = -k.s (s=dyT/dzT)
    ! Continent mask, multiply by dz/dy
    fluxm = dty / dtz * 500. * cuvsurcv * zmasqv
    ! slope limitation, multiply by kappa
    fluxm = -gmkappa * SIGN(MIN(ABS(fluxm), gm_smax * cv * cuvsurcv), dty)
    ! Zonal flux = 0. (temporary)
    fluxz(:) = 0.
    !  Vertical mass flux from mass budget (divergence of horiz fluxes)
    DO ij = iip2, ip1jm
      fluxv(ij) = fluxz(ij) - fluxz(ij - 1) - fluxm(ij) + fluxm(ij - iip1)
    ENDDO
    DO ij = iip1, ip1jmi1, iip1
      fluxv(ij + 1) = fluxv(ij + iip1)
    END DO
    !  vertical mass flux at Poles
    fluxv(1) = -SUM(fluxm(1:iim))
    fluxv(ip1jmp1) = SUM(fluxm(ip1jm - iim:ip1jm - 1))
    fluxv = fluxv

    ! Meridional heat fluxes
    DO ij = 1, ip1jm
      ! centered scheme
      fluxtm(ij, 1) = fluxm(ij) * (ts(ij + iip1, 1) + ts(ij, 1)) / 2.
      fluxtm(ij, 2) = -fluxm(ij) * (ts(ij + iip1, 2) + ts(ij, 2)) / 2.
    END DO

    ! Zonal heat fluxes
    ! Schema upstream
    DO ij = iip2, ip1jm
      IF (fluxz(ij)>=0.) THEN
        fluxtz(ij, 1) = fluxz(ij) * ts(ij, 1)
        fluxtz(ij, 2) = -fluxz(ij) * ts(ij + 1, 2)
      ELSE
        fluxtz(ij, 1) = fluxz(ij) * ts(ij + 1, 1)
        fluxtz(ij, 2) = -fluxz(ij) * ts(ij, 2)
      ENDIF
    ENDDO
    DO ij = iip1 * 2, ip1jmp1, iip1
      fluxtz(ij, :) = fluxtz(ij - iim, :)
    END DO

    ! Temperature tendency :
    DO ij = iip2, ip1jm
      ! divergence of horizontal heat fluxes
      dt(ij, :) = fluxtz(ij - 1, :) - fluxtz(ij, :) &
              + fluxtm(ij, :) - fluxtm(ij - iip1, :)
      ! + vertical heat flux (mass flux * T, upstream scheme)
      IF (fluxv(ij)>0.) THEN
        dt(ij, 1) = dt(ij, 1) + fluxv(ij) * ts(ij, 2)
        dt(ij, 2) = dt(ij, 2) - fluxv(ij) * ts(ij, 2)
      ELSE
        dt(ij, 1) = dt(ij, 1) + fluxv(ij) * ts(ij, 1)
        dt(ij, 2) = dt(ij, 2) - fluxv(ij) * ts(ij, 1)
      ENDIF
      ! divide by cell area
      dt(ij, :) = dt(ij, :) / aire(ij)
    END DO
    DO ij = iip1, ip1jmi1, iip1
      dt(ij + 1, :) = dt(ij + iip1, :)
    END DO
    ! Poles
    dt(1, :) = SUM(fluxtm(1:iim, :), dim = 1)
    IF (fluxv(1)>0.) THEN
      dt(1, 1) = dt(1, 1) + fluxv(1) * ts(1, 2)
      dt(1, 2) = dt(1, 2) - fluxv(1) * ts(1, 2)
    ELSE
      dt(1, 1) = dt(1, 1) + fluxv(1) * ts(1, 1)
      dt(1, 2) = dt(1, 2) - fluxv(1) * ts(1, 1)
    ENDIF
    dt(1, :) = dt(1, :) / apoln
    dt(ip1jmp1, :) = -SUM(fluxtm(ip1jm - iim:ip1jm - 1, :), dim = 1)
    IF (fluxv(ip1jmp1)>0.) THEN
      dt(ip1jmp1, 1) = dt(ip1jmp1, 1) + fluxv(ip1jmp1) * ts(ip1jmp1, 2)
      dt(ip1jmp1, 2) = dt(ip1jmp1, 2) - fluxv(ip1jmp1) * ts(ip1jmp1, 2)
    ELSE
      dt(ip1jmp1, 1) = dt(ip1jmp1, 1) + fluxv(ip1jmp1) * ts(ip1jmp1, 1)
      dt(ip1jmp1, 2) = dt(ip1jmp1, 2) - fluxv(ip1jmp1) * ts(ip1jmp1, 1)
    ENDIF
    dt(ip1jmp1, :) = dt(ip1jmp1, :) / apols
    dt(2:iip1, 1) = dt(1, 1)
    dt(2:iip1, 2) = dt(1, 2)
    dt(ip1jm + 1:ip1jmp1, 1) = dt(ip1jmp1, 1)
    dt(ip1jm + 1:ip1jmp1, 2) = dt(ip1jmp1, 2)

    ! T tendency back to 1D grid...
    CALL gr_dyn_fi(2, iip1, jjp1, dt, dt_phy)

  END SUBROUTINE slab_gmdiff

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE gr_fi_dyn(nfield, im, jm, pfi, pdyn)
    ! Transfer a variable from 1D "physics" grid to 2D "dynamics" grid
    USE lmdz_ssum_scopy, ONLY: scopy

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: im, jm, nfield
    REAL, INTENT(IN) :: pfi(klon_glo, nfield) ! on 1D grid
    REAL, INTENT(OUT) :: pdyn(im, jm, nfield) ! on 2D grid

    INTEGER :: i, j, ifield, ig

    DO ifield = 1, nfield
      ! Handle poles
      DO i = 1, im
        pdyn(i, 1, ifield) = pfi(1, ifield)
        pdyn(i, jm, ifield) = pfi(klon_glo, ifield)
      ENDDO
      ! Other points
      DO j = 2, jm - 1
        ig = 2 + (j - 2) * (im - 1)
        CALL SCOPY(im - 1, pfi(ig, ifield), 1, pdyn(1, j, ifield), 1)
        pdyn(im, j, ifield) = pdyn(1, j, ifield)
      ENDDO
    ENDDO ! of DO ifield=1,nfield

  END SUBROUTINE gr_fi_dyn

  SUBROUTINE gr_dyn_fi(nfield, im, jm, pdyn, pfi)
    ! Transfer a variable from 2D "dynamics" grid to 1D "physics" grid
    USE lmdz_ssum_scopy, ONLY: scopy
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: im, jm, nfield
    REAL, INTENT(IN) :: pdyn(im, jm, nfield) ! on 2D grid
    REAL, INTENT(OUT) :: pfi(klon_glo, nfield) ! on 1D grid

    INTEGER j, ifield, ig

    CHARACTER (len = 20) :: modname = 'slab_heat_transp'
    CHARACTER (len = 80) :: abort_message

    ! Sanity check:
    IF(klon_glo/=2 + (jm - 2) * (im - 1)) THEN
      abort_message = "gr_dyn_fi error, wrong sizes"
      CALL abort_physic(modname, abort_message, 1)
    ENDIF

    ! Handle poles
    CALL SCOPY(nfield, pdyn, im * jm, pfi, klon_glo)
    CALL SCOPY(nfield, pdyn(1, jm, 1), im * jm, pfi(klon_glo, 1), klon_glo)
    ! Other points
    DO ifield = 1, nfield
      DO j = 2, jm - 1
        ig = 2 + (j - 2) * (im - 1)
        CALL SCOPY(im - 1, pdyn(1, j, ifield), 1, pfi(ig, ifield), 1)
      ENDDO
    ENDDO

  END SUBROUTINE gr_dyn_fi

  SUBROUTINE  grad(klevel, pg, pgx, pgy)
    ! compute the covariant components pgx,pgy of the gradient of pg
    ! pgx = d(pg)/dx * delta(x) = delta(pg)
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: klevel
    REAL, INTENT(IN) :: pg((nbp_lon + 1) * nbp_lat, klevel)
    REAL, INTENT(OUT) :: pgx((nbp_lon + 1) * nbp_lat, klevel)
    REAL, INTENT(OUT) :: pgy((nbp_lon + 1) * (nbp_lat - 1), klevel)

    INTEGER :: l, ij
    INTEGER :: iim, iip1, ip1jm, ip1jmp1

    iim = nbp_lon
    iip1 = nbp_lon + 1
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1

    DO l = 1, klevel
      DO ij = 1, ip1jmp1 - 1
        pgx(ij, l) = pg(ij + 1, l) - pg(ij, l)
      ENDDO
      ! correction for pgx(ip1,j,l) ...
      ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
      DO ij = iip1, ip1jmp1, iip1
        pgx(ij, l) = pgx(ij - iim, l)
      ENDDO
      DO ij = 1, ip1jm
        pgy(ij, l) = pg(ij, l) - pg(ij + iip1, l)
      ENDDO
    ENDDO

  END SUBROUTINE grad

  SUBROUTINE diverg(klevel, x, y, div)
    ! computes the divergence of a vector field of components
    ! x,y. x and y being covariant components
    USE lmdz_ssum_scopy, ONLY: ssum

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: klevel
    REAL, INTENT(IN) :: x((nbp_lon + 1) * nbp_lat, klevel)
    REAL, INTENT(IN) :: y((nbp_lon + 1) * (nbp_lat - 1), klevel)
    REAL, INTENT(OUT) :: div((nbp_lon + 1) * nbp_lat, klevel)

    INTEGER :: l, ij
    INTEGER :: iim, iip1, iip2, ip1jm, ip1jmp1, ip1jmi1

    REAL :: aiy1(nbp_lon + 1), aiy2(nbp_lon + 1)
    REAL :: sumypn, sumyps

    iim = nbp_lon
    iip1 = nbp_lon + 1
    iip2 = nbp_lon + 2
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1
    ip1jmi1 = (nbp_lon + 1) * (nbp_lat - 1) - (nbp_lon + 1) ! = ip1jm - iip1

    DO l = 1, klevel
      DO ij = iip2, ip1jm - 1
        div(ij + 1, l) = &
                cvusurcu(ij + 1) * x(ij + 1, l) - cvusurcu(ij) * x(ij, l) + &
                        cuvsurcv(ij - iim) * y(ij - iim, l) - cuvsurcv(ij + 1) * y(ij + 1, l)
      ENDDO
      ! correction for div(1,j,l) ...
      ! ... div(1,j,l)= div(iip1,j,l) ...
      DO ij = iip2, ip1jm, iip1
        div(ij, l) = div(ij + iim, l)
      ENDDO
      ! at the poles
      DO ij = 1, iim
        aiy1(ij) = cuvsurcv(ij) * y(ij, l)
        aiy2(ij) = cuvsurcv(ij + ip1jmi1) * y(ij + ip1jmi1, l)
      ENDDO
      sumypn = SSUM(iim, aiy1, 1) / apoln
      sumyps = SSUM(iim, aiy2, 1) / apols
      DO ij = 1, iip1
        div(ij, l) = -sumypn
        div(ij + ip1jm, l) = sumyps
      ENDDO
      ! End (poles)
    ENDDO ! of DO l=1,klevel

    !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
    DO l = 1, klevel
      DO ij = iip2, ip1jm
        div(ij, l) = div(ij, l) * unsaire(ij)
      ENDDO
    ENDDO

  END SUBROUTINE diverg

  SUBROUTINE gr_v_scal(nx, x_v, x_scal)
    ! convert values from v points to scalar points on C-grid
    ! used to  compute unsfu, unseu (u points, but depends only on latitude)
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: nx ! number of levels or fields
    REAL, INTENT(IN) :: x_v((nbp_lon + 1) * (nbp_lat - 1), nx)
    REAL, INTENT(OUT) :: x_scal((nbp_lon + 1) * nbp_lat, nx)

    INTEGER :: l, ij
    INTEGER :: iip1, iip2, ip1jm, ip1jmp1

    iip1 = nbp_lon + 1
    iip2 = nbp_lon + 2
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1

    DO l = 1, nx
      DO ij = iip2, ip1jm
        x_scal(ij, l) = &
                (airev(ij - iip1) * x_v(ij - iip1, l) + airev(ij) * x_v(ij, l)) &
                        / (airev(ij - iip1) + airev(ij))
      ENDDO
      DO ij = 1, iip1
        x_scal(ij, l) = 0.
      ENDDO
      DO ij = ip1jm + 1, ip1jmp1
        x_scal(ij, l) = 0.
      ENDDO
    ENDDO

  END SUBROUTINE gr_v_scal

  SUBROUTINE gr_scal_v(nx, x_scal, x_v)
    ! convert values from scalar points to v points on C-grid
    ! used to compute wind stress at V points
    IMPLICIT NONE

    INTEGER, INTENT(IN) :: nx ! number of levels or fields
    REAL, INTENT(OUT) :: x_v((nbp_lon + 1) * (nbp_lat - 1), nx)
    REAL, INTENT(IN) :: x_scal((nbp_lon + 1) * nbp_lat, nx)

    INTEGER :: l, ij
    INTEGER :: iip1, ip1jm

    iip1 = nbp_lon + 1
    ip1jm = (nbp_lon + 1) * (nbp_lat - 1) ! = iip1*jjm

    DO l = 1, nx
      DO ij = 1, ip1jm
        x_v(ij, l) = &
                (cu(ij) * cvusurcu(ij) * x_scal(ij, l) + &
                        cu(ij + iip1) * cvusurcu(ij + iip1) * x_scal(ij + iip1, l)) &
                        / (cu(ij) * cvusurcu(ij) + cu(ij + iip1) * cvusurcu(ij + iip1))
      ENDDO
    ENDDO

  END SUBROUTINE gr_scal_v

  SUBROUTINE gr_scal_u(nx, x_scal, x_u)
    ! convert values from scalar points to U points on C-grid
    ! used to compute wind stress at U points
    USE lmdz_ssum_scopy, ONLY: scopy

    IMPLICIT NONE

    INTEGER, INTENT(IN) :: nx
    REAL, INTENT(OUT) :: x_u((nbp_lon + 1) * nbp_lat, nx)
    REAL, INTENT(IN) :: x_scal((nbp_lon + 1) * nbp_lat, nx)

    INTEGER :: l, ij
    INTEGER :: iip1, jjp1, ip1jmp1

    iip1 = nbp_lon + 1
    jjp1 = nbp_lat
    ip1jmp1 = (nbp_lon + 1) * nbp_lat ! = iip1*jjp1

    DO l = 1, nx
      DO ij = 1, ip1jmp1 - 1
        x_u(ij, l) = &
                (aire(ij) * x_scal(ij, l) + aire(ij + 1) * x_scal(ij + 1, l)) &
                        / (aire(ij) + aire(ij + 1))
      ENDDO
    ENDDO

    CALL SCOPY(nx * jjp1, x_u(1, 1), iip1, x_u(iip1, 1), iip1)

  END SUBROUTINE gr_scal_u

END MODULE slab_heat_transp_mod