*IDENT wmc020801
*DECLARE ICEFREED
*/
*/ Problems with this: it should be in vn4.5 as a deck, but doesn't
*/ show up properly when linked
*/
*/ In the best of traditions, fudge this by adding it into ICEFREED
*/
*I ICEFREED.255
C ******************************COPYRIGHT******************************
C (c) CROWN COPYRIGHT 1998, METEOROLOGICAL OFFICE, All Rights Reserved.
C
C Use, duplication or disclosure of this code is subject to the
C restrictions as set forth in the contract.
C
C                Meteorological Office
C                London Road
C                BRACKNELL
C                Berkshire UK
C                RG12 2SZ
C
C If no contract has been raised with this copy of the code, the use,
C duplication or disclosure of it is strictly prohibited.  Permission
C to do so must first be obtained in writing from the Head of Numerical
C Modelling at the above address.
C ******************************COPYRIGHT******************************
!+Convergence stopping routine for sea-ice
!
! Subroutine Interface:
       subroutine cnvstop(
*CALL ARGOINDX
     & imt,imtm1,imtm2,jmt,jmtm1,jmtm2,
     & uice,
     & vice,
     & hice,
     & radius,
     & dphi,
     & dlambda,
     & phit,
     & hicemax,hicemid
     & )
      IMPLICIT NONE
!
! Description:
! This subroutine restricts the flow of sea-ice up a gradient in HICE
! (GBM sea-ice depth), as pseudo-rheological resistance to convergence.
! The flow parallel to grad(HICE) is impeded. Where HICE < HICEMID,
! the flow is uneffected, and where HICE > HICEMAX the flow is stopped
! completely, with ! a linear ramping of the change to the velocity
! between HICEMID and HICEMAX. It is set up to work with the freedrift
! scheme, being called in ICEFREEDR, but could equally well work with
! the simple ice advection scheme, replacing ORH3F403.76 to
! ICEDRIFT.271 (at vn4.4) of ICEDRIFT.
!
! Method:
! Looping over velocity points, the maximum of surrounding values of
! HICE is taken and labelled HICE_U. This value dictates whether, and
! how, the ice velocity will be restricted. If HICE_U > HICEMID the
! components of grad(HICE) and U.grad(HICE) are calculated. If
! U.grad(HICE) is positive (i.e.the ice is deep and converging) then
! the component of U parallel to grad(HICE) is scaled s.t. it is
! unchanged where HICE_U = HICEMID, and zeroed where HICE_U = HICEMAX,
! the scaling varying linearly between the two states.
!
! Current Code Owner: Doug Cresswell
!
! History:
! Version   Date     Comment
! -------   ----     -------
!   4.5   25.9.98    New deck. Doug Cresswell and Jonathan Gregory.
!
! Code Description:
!   Language: FORTRAN 77 + common extensions.
!   This code is written to UMDP3 v6 programming standards.
C
*CALL OTIMER
*CALL TYPOINDX
*CALL CNTLOCN
*CALL OARRYSIZ
       integer
C*CALL ARGSIZE
C      integer
     & imt            ! number of tracer columns.
     &,imtm1          ! number of tracer columns - 1
     &,imtm2          ! number of tracer columns - 2
     &,jmt            ! number of tracer rows.
     &,jmtm1          ! number of velocity rows.
     &,jmtm2          ! number of velocity rows - 1
      real
     & uice(imt,jmtm1)  ! inout   zonal sea ice velocity.
     &,vice(imt,jmtm1)  ! inout   meridional sea ice velocity.
     &,hice(imt,jmt)    ! in ice thickness.
     &,radius           ! in radius of the earth in metres.
     &,dphi(jmt)        ! in meridional grid spacing in radians.
     &,dlambda(imt)     ! in zonal grid spacing in radians.
     &,phit(jmt)        ! in    latitude of mass rows in radians.
     &,hmask(imt,jmt)   ! in 1.0 for ha land 0.0 for sea.
     &,umask(imt,jmtm1) ! in 1.0 for uv land 0.0 for sea.
     &,hicemax          ! in max hice at which convergence is allowed.
     &,hicemid          ! in min hice at which convergence is impeded.
C Local Variables
      integer
     & i,j                 !Looping variables
      real
     & facex,facey,faceyp  !Lengths of the E(&W),S and N sides
                           !of the gridbox.
     &,gradhu(imt,jmtm1)   ! U cpt of the grad h (ice depth)
     &,gradhv(imt,jmtm1)   ! V cpt of the grad h (ice depth)
     &,modgradhsq(imt,jmtm1)   ! Mod(Grad h), squared.
     &,udotgradh(imt,jmtm1)    ! Scalar product of u and grad h
     &,hice_u(imt,jmtm1)   ! h interpolated onto the U grid
     &,recip_hmax_m_hmid   !1/(hicemax-hicemin) (max=100m-1)
     &,b                   !Scaling factor for ice between hicemid
                           !and hicemax.
C Start executable code
      if (L_OTIMER) call timer('cnvstop',3)
*IF DEF,MPP
      CALL SWAPBOUNDS(UICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
      CALL SWAPBOUNDS(VICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
*/ MJR to JMT?
*/     CALL SWAPBOUNDS(HICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
      CALL SWAPBOUNDS(HICE,IMT,JMT,O_EW_HALO,O_NS_HALO,1)
*ENDIF
      recip_hmax_m_hmid=1.0/max((hicemax-hicemid),0.01)
C Start looping over u points
      do j=J_1,J_jmt
        do i=2,imtm1
C Initialise arrays to zero, to be partially overwritten
C (Note that these arrays are not made cyclic, as they are not output)
            gradhu(i,j)=0.0
            gradhv(i,j)=0.0
            modgradhsq(i,j)=0.0
            udotgradh(i,j)=0.0
C Take hice on the u grid to be the maximum of surrounding hice's. This
C prevents advection into deep areas up steep gradients in hice.
            hice_u(i,j)=max(hice(i,j),hice(i+1,j)
     &                      ,hice(i,j+1),hice(i+1,j+1))
C If ice is deep then procede.
            if(hice_u(i,j).gt.hicemid) then
C Calculate gridbox side lengths
              facex  = radius*dphi(j)
              facey  = (radius*cos(phit(j)))*dlambda(i+1)
              faceyp = (radius*cos(phit(j+1)))*dlambda(i+1)
C Calculate Grad h (u and v cpts)
              gradhu(i,j)=0.5*(hice(i+1,j+1)-hice(i,j+1))/faceyp
     &                   +0.5*(hice(i+1,j)-hice(i,j))/facey
              gradhv(i,j)=0.5*((hice(i,j+1)-hice(i,j))
     &                        +(hice(i+1,j+1)-hice(i+1,j)))/facex
C Calculate scalar product of U and Grad h
              udotgradh(i,j)=uice(i,j)*gradhu(i,j)
     &                       +vice(i,j)*gradhv(i,j)
C If U is acting up the ice depth (h) gradient then
C set the cpt of U parallel to Grad h to zero.
              if(udotgradh(i,j).gt.0.0) then
                  modgradhsq(i,j)=gradhu(i,j)*gradhu(i,j)
     &                           +gradhv(i,j)*gradhv(i,j)
                  if(hice_u(i,j).gt.hicemax)then
                    uice(i,j)=uice(i,j)-(udotgradh(i,j)*gradhu(i,j))
     &                                   /modgradhsq(i,j)
                    vice(i,j)=vice(i,j)-(udotgradh(i,j)*gradhv(i,j))
     &                                   /modgradhsq(i,j)
                  else if(hice_u(i,j).gt.hicemid)then
                    b=(hice_u(i,j)-hicemid)*recip_hmax_m_hmid
                    uice(i,j)=uice(i,j)-b*(udotgradh(i,j)*gradhu(i,j))
     &                                   /modgradhsq(i,j)
                    vice(i,j)=vice(i,j)-b*(udotgradh(i,j)*gradhv(i,j))
     &                                   /modgradhsq(i,j)
                  endif
              endif
          endif
C Looping - need to sort out cyclic points
        enddo
      enddo
       if(L_OCYCLIC) then
         do j=J_1,J_jmt
           uice(1,j)=uice(imtm1,j)
           uice(imt,j)=uice(2,j)
           vice(1,j)=vice(imtm1,j)
           vice(imt,j)=vice(2,j)
         enddo
       endif
      CALL SWAPBOUNDS(UICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
      CALL SWAPBOUNDS(VICE,IMT,JMTM1,O_EW_HALO,O_NS_HALO,1)
      if (L_OTIMER) call timer('cnvstop',4)
      end