*IDENT STERNc2d
*/
*/ Author: WMC
*/
*/ Purpose:       implements a parametrisation of open-water creation by
*/                shear etc, as given in Stern et al., JGR 100, C10, p20601, 1995
*/                note that (in accordance with std practice) the open-water
*/                created is /20.
*/ Category:      seaice
*/ Prerequisites: EVP seaice mod evp-new1.upd
*/ Revisions:     wmc 2003/03/03 
*/                - correct code (and diag) for when aice would become -ve
*/                - arbitrarily assert no ridging for aice < 0.5
*/                  (this makes error condition very very unlikely)
*/                - make default that 1/2 energy goes into ridging
*/
*/ Stern eq 16 says: phi_thin = e_1_dot + 1/2*(delta - e_1_dot)*W
*/ Since our h_0=0 (ie, thin ice is open water), W==1
*/
*/ Now, this term *replaces* the usual term e_1_dot, and I *think* that
*/ e_1_dot is included, by default, as the divergence. So I *think* that
*/ the term we wish to add is:
*/   
*/   phi_t = 1/2*(delta - e_1_dot)*W
*/
*/ Now, delta is available from EVP as deltan/e/w/s (each grid square is
*/ in triangles) which I presume is best averaged. Also e_1_dot is
*/ available as xi11n/e/s/w, ditto.
*/
*/ Mod evp-new1.upd arranges for this variable to be created and
*/ passed up into ice_ctl. This mod uses it.
*/
*/ Note: phi_t requires mult by the timestep to be used. I have hardcoded
*/ 1 hour (3600) in this mod. Sorry.
*/
*/ Note: there may well be a case for using the alternative form by Hibler,
*/ (see-also CICE code) which allows you to use the divergence as actually
*/ implemented in the advection routine rather than its theoretical value
*/ (Lipscomb says: no: don't do this: and has reverted from doing so).

*DECLARE ICE_CTL

*/ Now use phi_t. This is analagous to the old c2d mod

*I ORH1F305.5732
C
C WMC. improved conc -> depth parametrisation.
C
C      Since icemass = dens * htrue * aice = dens * hice,
C      (and obviously we want to conserve mass),
C      all we have to do is adjust aice downwards.
C

C First, special code to record where errors occurred - 
C errors being where we're trying to make the conc lt 0
C Should really be integer but...
C Revision: where aice is zero, don't count
       IF (SF(163,32)) THEN
         DO J = 1, JMT
          DO i=1, IMT
             aice_ctod(i,j)=phi_t(i,j)*3600/20.
             IF (aice(i,j) .gt. 0.5 .and.
     &           aice(i,j) .lt. aice_ctod(i,j))
     &       THEN
               aice_ctod_err(i,j)=1.0
             ELSE
               aice_ctod_err(i,j)=0.0
             ENDIF
          END DO
         END DO
         CALL COPYODIAGL(IMT,JMT,1,.TRUE.,0.
     &         ,aice_ctod_err,ocean,STASHWORK(SI(163,32,im_index)))
       ENDIF

C Now actually do the calculation
       DO J = 1, JMT
         DO i=1, IMT

C Note that we have naughtily hardcoded the timestep in seconds
             aice_ctod(i,j)=phi_t(i,j)*3600.0/20.
C Could this make aice -ve? Lets test...
             IF (aice(i,j) .gt. 0.5 .and.
     &           aice(i,j) .lt. aice_ctod(i,j))
     &       THEN
	       print *,'Oops: -ve aice in c2d: ',
     &                 i,j,aice(i,j),aice_ctod(i,j)
             ELSE
               aice(i,j)=aice(i,j)-aice_ctod(i,j)
             ENDIF

         END DO
       END DO
       IF (SF(161,32)) CALL COPYODIAGL(IMT,JMT,1,.TRUE.,0.
     &              ,aice_ctod,ocean,STASHWORK(SI(161,32,im_index)))

*B OJG3F403.13  
C WMC. For stash items 32.161/3 (evp_vn3)
C      Could dimension as imt_ice, jmt_ice but why bother? You'll only ask
C      for this with seaice running...
      REAL                           
     & aice_ctod(imt,jmt),                         ! Fraction lost by c2d
     & aice_ctod_err(imt,jmt)                      ! Set to 1 when fraction would be < 0