!        Generated by TAPENADE     (INRIA, Tropics team)
!  Version 2.1 - (Id: 1.33 vmp Stable - Wed Jan 12 13:35:07 MET 2005)
!  
! $Id: ice_ocean.F,v 1.9 2003/11/22 21:34:50 lipscomb Exp $
!=======================================================================
!BOP
!
! !MODULE: ice_ocean - ocean mixed layer internal to sea ice model
!
! !DESCRIPTION:
!
! Ocean mixed layer calculation (internal to sea ice model).
! Allows heat storage in ocean for uncoupled runs.
!
! !REVISION HISTORY:
!
! authors:   John Weatherly, CRREL
!            C.M. Bitz, UW
!            Elizabeth C. Hunke, LANL
!            Bruce P. Briegleb, NCAR
!            William H. Lipscomb, LANL
!
! !INTERFACE:
!
MODULE ICE_OCEAN_DV
  USE ice_albedo_dv
  USE ice_atmo_dv
  USE ice_calendar_dv
  USE ice_constants_dv
  USE ice_domain_dv
  USE ice_flux_dv
  USE ice_grid_dv
  USE ice_kinds_mod_dv
  USE ice_model_size_dv
  USE ice_state_dv
  REAL(KIND=DBL_KIND) :: cphm, cphmd(nbdirsmax), hmix, hmixd(nbdirsmax)
  LOGICAL(KIND=LOG_KIND) :: oceanmixed_ice

CONTAINS
!  Differentiation of mixed_layer in forward (tangent) mode: (multi-directional mode)
!   variations  of output variables: lhcoef frzmlt shcoef sst
!   with respect to input variables: cphm aice qa lhcoef frzmlt
!                swvdf wind swidf rhoa swvdr swidr shcoef pott
!                flw sst tf emissivity
  SUBROUTINE MIXED_LAYER_DV(nbdirs)
    IMPLICIT NONE
    INCLUDE 'DIFFSIZES.inc'
!  Hint: nbdirsmax should be the maximum number of differentiation directions
    INTEGER :: nbdirs
    REAL(KIND=DBL_KIND) :: delq(ilo:ihi, jlo:jhi), delqd(nbdirsmax, ilo:&
&    ihi, jlo:jhi), delt(ilo:ihi, jlo:jhi), deltd(nbdirsmax, ilo:ihi, &
&    jlo:jhi), dummy1(ilo:ihi, jlo:jhi), dummy1d(nbdirsmax, ilo:ihi, jlo&
&    :jhi), dummy2(ilo:ihi, jlo:jhi), dummy2d(nbdirsmax, ilo:ihi, jlo:&
&    jhi), dummy3(ilo:ihi, jlo:jhi), dummy3d(nbdirsmax, ilo:ihi, jlo:jhi&
&    ), dummy4(ilo:ihi, jlo:jhi)
    REAL(KIND=DBL_KIND) :: flh, flhd(nbdirsmax), flwup, flwupd(nbdirsmax&
&    ), fsh, fshd(nbdirsmax), ft, qdp, swabs, swabsd(nbdirsmax), tsfk, &
&    tsfkd(nbdirsmax)
    INTEGER(KIND=INT_KIND) :: i, j
    INTEGER :: nd
    REAL :: x1, x1d(nbdirsmax)
    INTRINSIC MAX, MIN
    DO nd=1,nbdirs
      vatmd(nd, :, :) = 0.0
      delqd(nd, :, :) = 0.0
      deltd(nd, :, :) = 0.0
      dummy3d(nd, :, :) = 0.0
      dummy2d(nd, :, :) = 0.0
      dummy1d(nd, :, :) = 0.0
      uatmd(nd, :, :) = 0.0
    END DO
!
! !USES:
! 
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! horizontal indices
!
! potential temperature difference   (K)
! specific humidity difference   (kg/kg)
! dummy arrays
!
! surface temperature (K)
! sensible heat flux  (W/m^2)
! latent heat flux    (W/m^2)
! surface absorbed shortwave heat flux (W/m^2)
! long-wave upward heat flux  (W/m^2)
! deep ocean heat flux
! fraction reduction of positive qdp
!
    CALL ATMO_BOUNDARY_LAYER_DV(1, 'ocn', sst, sstd, dummy1, dummy1d, &
&                          dummy2, dummy2d, dummy3, dummy3d, dummy4, &
&                          delt, deltd, delq, delqd, nbdirs)
!
    DO j=jlo,jhi
      DO i=ilo,ihi
! tmask
        IF (tmask(i, j)) THEN
! specify deep ocean heat flux as constant for now
!
! negative upward
! ocean surface temperature in Kelvin
          qdp = -c2
!
! shortwave radiative flux
          tsfk = sst(i, j) + tffresh
!
! longwave radiative flux
          swabs = (c1-albocn)*(swvdr(i, j)+swvdf(i, j)+swidr(i, j)+swidf&
&            (i, j))
!
! downward latent and sensible heat fluxes
          flwup = -(emissivity*stefan_boltzmann*tsfk**4)
!
          flh = lhcoef(i, j)*delq(i, j)
! first, compute sst change due to exchange with atm/ice above
          fsh = shcoef(i, j)*delt(i, j)
!
! computed T change due to exchange with deep layers:
          sst(i, j) = sst(i, j) + (fsh+flh+flwup+flw(i, j)+swabs)*(c1-&
&            aice(i, j))*dt/cphm
!
! compute potential to freeze or melt ice
          sst(i, j) = sst(i, j) - qdp*dt/cphm
          DO nd=1,nbdirs
            swabsd(nd) = (c1-albocn)*(swvdrd(nd, i, j)+swvdfd(nd, i, j)+&
&              swidrd(nd, i, j)+swidfd(nd, i, j))
            tsfkd(nd) = sstd(nd, i, j)
            fshd(nd) = shcoefd(nd, i, j)*delt(i, j) + shcoef(i, j)*deltd&
&              (nd, i, j)
            flhd(nd) = lhcoefd(nd, i, j)*delq(i, j) + lhcoef(i, j)*delqd&
&              (nd, i, j)
            flwupd(nd) = -(stefan_boltzmann*(emissivityd(nd)*tsfk**4+&
&              emissivity*4*tsfk**3*tsfkd(nd)))
            sstd(nd, i, j) = sstd(nd, i, j) + (dt*((fshd(nd)+flhd(nd)+&
&              flwupd(nd)+flwd(nd, i, j)+swabsd(nd))*(c1-aice(i, j))-(&
&              fsh+flh+flwup+flw(i, j)+swabs)*aiced(nd, i, j))*cphm-(fsh&
&              +flh+flwup+flw(i, j)+swabs)*(c1-aice(i, j))*dt*cphmd(nd))&
&              /cphm**2
            sstd(nd, i, j) = sstd(nd, i, j) + qdp*dt*cphmd(nd)/cphm**2
            frzmltd(nd, i, j) = ((tfd(nd, i, j)-sstd(nd, i, j))*cphm+(tf&
&              (i, j)-sst(i, j))*cphmd(nd))/dt
          END DO
!
          frzmlt(i, j) = (tf(i, j)-sst(i, j))*cphm/dt
          IF (frzmlt(i, j) .LT. -c1000) THEN
            x1 = -c1000
            DO nd=1,nbdirs
              x1d(nd) = 0.0
            END DO
          ELSE
            DO nd=1,nbdirs
              x1d(nd) = frzmltd(nd, i, j)
            END DO
            x1 = frzmlt(i, j)
          END IF
          IF (x1 .GT. c1000) THEN
            DO nd=1,nbdirs
              frzmltd(nd, i, j) = 0.0
            END DO
            frzmlt(i, j) = c1000
          ELSE
            DO nd=1,nbdirs
              frzmltd(nd, i, j) = x1d(nd)
            END DO
            frzmlt(i, j) = x1
          END IF
!
          IF (sst(i, j) .LE. tf(i, j)) THEN
            DO nd=1,nbdirs
              sstd(nd, i, j) = tfd(nd, i, j)
            END DO
            sst(i, j) = tf(i, j)
          END IF
        END IF
      END DO
    END DO
  END SUBROUTINE MIXED_LAYER_DV
  SUBROUTINE MIXED_LAYER()
    IMPLICIT NONE
    REAL(KIND=DBL_KIND) :: delq(ilo:ihi, jlo:jhi), delt(ilo:ihi, jlo:jhi&
&    ), dummy1(ilo:ihi, jlo:jhi), dummy2(ilo:ihi, jlo:jhi), dummy3(ilo:&
&    ihi, jlo:jhi), dummy4(ilo:ihi, jlo:jhi)
    REAL(KIND=DBL_KIND) :: flh, flwup, fsh, ft, qdp, swabs, tsfk
    INTEGER(KIND=INT_KIND) :: i, j
    REAL :: x1
    INTRINSIC MAX, MIN
!
! !USES:
! 
!
!
! !INPUT/OUTPUT PARAMETERS:
!
!EOP
!
! horizontal indices
!
! potential temperature difference   (K)
! specific humidity difference   (kg/kg)
! dummy arrays
!
! surface temperature (K)
! sensible heat flux  (W/m^2)
! latent heat flux    (W/m^2)
! surface absorbed shortwave heat flux (W/m^2)
! long-wave upward heat flux  (W/m^2)
! deep ocean heat flux
! fraction reduction of positive qdp
!
    CALL ATMO_BOUNDARY_LAYER(1, 'ocn', sst, dummy1, dummy2, dummy3, &
&                       dummy4, delt, delq)
!
    DO j=jlo,jhi
      DO i=ilo,ihi
! tmask
        IF (tmask(i, j)) THEN
! specify deep ocean heat flux as constant for now
!
! negative upward
! ocean surface temperature in Kelvin
          qdp = -c2
!
! shortwave radiative flux
          tsfk = sst(i, j) + tffresh
!
! longwave radiative flux
          swabs = (c1-albocn)*(swvdr(i, j)+swvdf(i, j)+swidr(i, j)+swidf&
&            (i, j))
!
! downward latent and sensible heat fluxes
          flwup = -(emissivity*stefan_boltzmann*tsfk**4)
!
          flh = lhcoef(i, j)*delq(i, j)
! first, compute sst change due to exchange with atm/ice above
          fsh = shcoef(i, j)*delt(i, j)
!
! computed T change due to exchange with deep layers:
          sst(i, j) = sst(i, j) + (fsh+flh+flwup+flw(i, j)+swabs)*(c1-&
&            aice(i, j))*dt/cphm
!
! compute potential to freeze or melt ice
          sst(i, j) = sst(i, j) - qdp*dt/cphm
!
          frzmlt(i, j) = (tf(i, j)-sst(i, j))*cphm/dt
          IF (frzmlt(i, j) .LT. -c1000) THEN
            x1 = -c1000
          ELSE
            x1 = frzmlt(i, j)
          END IF
          IF (x1 .GT. c1000) THEN
            frzmlt(i, j) = c1000
          ELSE
            frzmlt(i, j) = x1
          END IF
!
          IF (sst(i, j) .LE. tf(i, j)) sst(i, j) = tf(i, j)
        END IF
      END DO
    END DO
  END SUBROUTINE MIXED_LAYER
END MODULE ICE_OCEAN_DV