--- MITgcm/pkg/seaice/seaice_jacvec.F	2012/11/06 13:09:30	1.2
+++ MITgcm/pkg/seaice/seaice_jacvec.F	2017/05/23 16:24:46	1.8
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.2 2012/11/06 13:09:30 mlosch Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.8 2017/05/23 16:24:46 mlosch Exp $
 C $Name:  $
 
 #include "SEAICE_OPTIONS.h"
@@ -56,9 +56,7 @@
       _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 
-#if ( defined (SEAICE_CGRID) && \
-      defined (SEAICE_ALLOW_JFNK) && \
-      defined (SEAICE_ALLOW_DYNAMICS) )
+#ifdef SEAICE_ALLOW_JFNK
 C     Local variables:
       _RL utp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL vtp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
@@ -71,19 +69,50 @@
       _RL epsilon, reps
 CEOP
 C     Instructions for using TAF or TAMC to generate exact Jacobian times
-C     vector operations:
+C     vector operations (if SEAICE_ALLOW_MOM_ADVECTION is defined, the
+C     file list also needs to include seaice_mom_advection.f, 
+C     mom_calc_hfacz.f, mom_calc_ke.f, mom_calc_relvort3.f, 
+C     mom_vi_u_coriolis.f, mom_vi_u_coriolis_c4.f, mom_vi_u_grad_ke.f,
+C     mom_vi_v_coriolis.f, mom_vi_v_coriolis_c4.f, mom_vi_v_grad_ke.f
+C     plus flow information for diagnostics_fill.f:
+CCCCCCCADJ SUBROUTINE DIAGNOSTICS_FILL INPUT  = 1,2,3,4,5,6,7,8
+CCCCCCCADJ SUBROUTINE DIAGNOSTICS_FILL OUTPUT =
+C     )
 C
 C     1. make small_f
-C     2. cat seaice_calc_residual.f seaice_oceandrag_coeffs.f seaice_calc_strainrates.f seaice_calc_viscosities.f seaice_calc_rhs.f seaice_calc_lhs.f > taf_input.f
-C     3. staf -v1 -forward -toplevel seaice_calc_residual -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
+C     2. cat seaice_calc_residual.f seaice_oceandrag_coeffs.f \
+C        seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f \
+C        seaice_calc_strainrates.f seaice_calc_viscosities.f \
+C        seaice_calc_rhs.f seaice_calc_lhs.f > taf_input.f
+C     3. staf -v1 -forward -toplevel seaice_calc_residual \
+C             -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
 C     4. insert content of taf_input_ftl.f at the end of this file
 C     5. add the following code and comment out the finite difference code
+C
+C     Instruction for using TAF 2.4 and higher (or staf with default -v2
+C     starting with version 2.0):
+C
+C     1. make small_f
+C     2. files="seaice_calc_residual.f seaice_oceandrag_coeffs.f \
+C               seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f \
+C               seaice_calc_strainrates.f seaice_calc_viscosities.f \
+C               seaice_calc_rhs.f seaice_calc_lhs.f"
+C     3. staf -forward -toplevel seaice_calc_residual \
+C             -input uIceLoc,viceLoc -output uIceRes,vIceRes $files
+C     4. copy files seaice_*_tl.f to the corresponding seaice_*.f files, 
+C        e.g. with this bash script:
+C     for file in $files; do 
+C       nfile=`echo $file | awk -F. '{printf "%s_tl.f", $1}'`; 
+C       \cp -f $nfile $file
+C     done
+C     5. add the following code, change "call g_seaice_calc_residual" 
+C        to "call seaice_calc_residual_tl", and comment out the finite 
+C        difference code
 CML      _RL g_duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 CML      _RL g_dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 CML      _RL g_uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 CML      _RL g_vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
 CML
-CML      IF (.false.) then
 CMLC     Initialise
 CML      DO bj=myByLo(myThid),myByHi(myThid)
 CML       DO bi=myBxLo(myThid),myBxHi(myThid)
@@ -103,6 +132,10 @@
 CML      CALL G_SEAICE_CALC_RESIDUAL( uIce, g_duice, vIce, 
 CML     $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter, 
 CML     $kryloviter, mytime, myiter, mythid )
+CMLCML      For staf -v2 replace the above with the below call 
+CMLCML      CALL SEAICE_CALC_RESIDUAL_TL( uIce, g_duice, vIce, 
+CMLCML     $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter, 
+CMLCML     $kryloviter, mytime, myiter, mythid )
 CML
 CML      DO bj=myByLo(myThid),myByHi(myThid)
 CML       DO bi=myBxLo(myThid),myBxHi(myThid)
@@ -116,7 +149,7 @@
 CML      ENDDO
 
 C     Initialise
-      epsilon = 1. _d -06
+      epsilon = SEAICE_JFNKepsilon
       reps    = 1. _d 0/epsilon
 
       DO bj=myByLo(myThid),myByHi(myThid)
@@ -151,7 +184,7 @@
        ENDDO
       ENDDO
 
-#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */
+#endif /* SEAICE_ALLOW_JFNK */
 
       RETURN
       END