C $Header: /u/gcmpack/MITgcm/pkg/streamice/streamice_invert_surf_forthick.F,v 1.4 2016/11/29 12:38:43 dgoldberg Exp $
C $Name:  $

#include "CPP_OPTIONS.h"
#include "STREAMICE_OPTIONS.h"

      SUBROUTINE STREAMICE_INVERT_SURF_FORTHICK (
     O            H,
     I            S,
     I            R,
     I            delta,
     I            myThid)

!      This S/R finds a thickness (H) that gives surf elev S with bed R


#include "SIZE.h"
#include "GRID.h"
#include "SET_GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "STREAMICE.h"

#ifdef ALLOW_OPENAD
      use OAD_tape
      use OAD_rev
      use OAD_cp
#endif

      _RL H(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL S(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL R(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL DELTA
      INTEGER myThid
#ifdef ALLOW_OPENAD
      type(active) :: ETA_GL_STREAMICE
      type(active) :: X,Y0
      type(modeType) :: our_orig_mode
#endif
      

#ifdef ALLOW_STREAMICE

      _RL resid, f, fp, hf, htmp
      INTEGER  i, j, bi, bj, ITER
      _RL ETA_GL_PRIME_STREAMICE
#ifndef ALLOW_OPENAD
      _RL ETA_GL_STREAMICE
      EXTERNAL 
#endif
!      EXTERNAL PHI_GL_STREAMICE_PRIME

        DO bj=myByLo(myThid), myByHi(myThid)
         DO bi=myBxLo(myThid), myBxHi(myThid)
          DO j = 1,sNy
           DO i = 1,sNx
            IF (STREAMICE_hmask(i,j,bi,bj).eq.1.0) THEN
             hf = (-1. _d 0) * R(i,j,bi,bj) /
     &           (1. _d 0 - delta)

             IF (S(i,j,bi,bj) .gt. delta*HF) THEN
              htmp = S(i,j,bi,bj)-R(i,j,bi,bj)
             ELSE
              htmp = S(i,j,bi,bj)/delta
             ENDIF

             IF (streamice_smooth_gl_width.gt.0.) THEN

              RESID=1. _d 0

              DO ITER=1,20
               IF ((RESID .gt. .005) .and.
     &           ( STREAMICE_hmask(i,j,bi,bj).eq.1.0)) THEN

                hf = (-1. _d 0) * R(i,j,bi,bj) /
     &           (1. _d 0 - delta)

                IF (S(i,j,bi,bj) .gt. delta*HF) THEN
                 htmp = S(i,j,bi,bj)-R(i,j,bi,bj)
                ELSE
                 htmp = S(i,j,bi,bj)/delta
                ENDIF

#ifdef ALLOW_OPENAD

                our_orig_mode = our_rev_mode
                our_rev_modearg_store=.FALSE.
                our_rev_modearg_restore=.FALSE.
                our_rev_modeplain=.TRUE.
                our_rev_modetape=.FALSE.
                our_rev_modeadjoint=.FALSE.

                Xv = htmp-HF
                Y0v = delta*HF

                CALL OPENAD_OAD_S_ETA_GL_STREAMICE(
     &             X,
     &             delta,
     &             1. _d 0,
     &             Y0,
     &             streamice_smooth_gl_width,
     &             ETA_GL_STREAMICE)

                RESID = ETA_GL_STREAMICEv
                our_rev_mode = our_orig_mode

#else
                RESID = ETA_GL_STREAMICE (
     &             htmp-HF,
     &             delta,
     &             1. _d 0,
     &             delta*HF,
     &             streamice_smooth_gl_width)
#endif
                RESID = RESID - S(i,j,bi,bj)
                FP = ETA_GL_PRIME_STREAMICE (
     &             htmp-HF,
     &             delta,
     &             1. _d 0,
     &             delta*HF,
     &             streamice_smooth_gl_width)
                Htmp = HTMP-RESID/FP
               ENDIF
              ENDDO
             ENDIF
             H(i,j,bi,bj) = Htmp
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO


#endif
      RETURN
      END