CALVESSEL.f

      SUBROUTINE calvessel(ICEG)
! Algorithm for calculating ice growth rate
      use vrbls2d, only: sst, u10h, v10h, tshltr
      use masks, only: sm, sice
      use ctlblk_mod, only: jsta, jend, im, spval
!-------------------------------------------
      implicit none
      integer i, j
      real tsfc_c,tshltr_c,sst_c
      real, parameter :: c2k=273.15
      real, dimension(im,jsta:jend) :: pr, spd10
      real,intent(out) ::  iceg(im,jsta:jend)

!      allocate (thsfc(im,jsta:jend),tsfc(im,jsta:jend))

      DO j=jsta,jend
        DO i=1,im
!   CALCULATE SPEED
          spd10(i,j)=sqrt(u10h(i,j)**2+v10h(i,j)**2)
          if (spd10(i,j)>50) then
            iceg(i,j)=0.
            cycle
          endif

! Reverse of land mask use le instead of ge from original code
!!  MASK CHECK
          if((sice(i,j)>=0.5).or.(sm(i,j)<=0.5)) then
            iceg(i,j)=0.
            cycle
          endif

!!! CHANGE TEMP to FROM K to C
!!! TEMPERATURE CHECK
          sst_c=sst(i,j)-c2k  
          tshltr_c=tshltr(i,j)-c2k
          if((sst_c<-1.7).OR. &
             (sst_c>12.0)) then
            iceg(i,j)=0.
            cycle
          endif

          if((tshltr_c>0.).OR. &
             (tshltr_c<-40.)) then
            iceg(i,j)=0.
            cycle
          endif

!  CALCULATE ICE GROWTH
          pr(i,j)=spd10(i,j)*(-1.7-tshltr_c)/(1.+.4*(sst_c+1.7))
          iceg(i,j)=(2.73e-02)*pr(i,j)+(2.91e-04)*pr(i,j)*pr(i,j) &
                   +(1.84e-06)*pr(i,j)**3

!! ICE GROWTH CHECK
          if (iceg(i,j)<0.) THEN
            iceg(i,j)=0.
          else
! Convert to m/s from cm/hr
            iceg(i,j)=(1./3.6e+05)*iceg(i,j)
          endif

        ENDDO
      ENDDO

      END