pkg/seaice/seaice_jacvec.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.3 2012/11/06 13:10:14 mlosch Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CBOP
C     !ROUTINE: SEAICE_JACVEC
C     !INTERFACE:
      SUBROUTINE SEAICE_JACVEC( 
     I     uIceLoc, vIceLoc, uIceRes, vIceRes,
     U     duIce, dvIce,  
     I     newtonIter, krylovIter, myTime, myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SEAICE_JACVEC
C     | o For Jacobian-free Newton-Krylov solver compute
C     |   Jacobian times vector by finite difference approximation
C     *==========================================================*
C     | written by Martin Losch, Oct 2012
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: my Thread Id. number
C     newtonIter :: current iterate of Newton iteration
C     krylovIter :: current iterate of Krylov iteration
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
      INTEGER newtonIter
      INTEGER krylovIter
C     u/vIceLoc :: local copies of the current ice velocity
      _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     u/vIceRes :: initial residual of this Newton iterate
      _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
C     du/vIce   :: correction of ice velocities
      _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) )
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)
C     u/vIceResP :: residual computed with u/vtp
      _RL uIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

C     i,j,bi,bj :: loop indices
      INTEGER i,j,bi,bj
      _RL epsilon, reps
CEOP
C     Instructions for using TAF or TAMC to generate exact Jacobian times
C     vector operations:
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     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. staf -forward -toplevel seaice_calc_residual -input uIceLoc,viceLoc -output uIceRes,vIceRes seaice_calc_residual.f seaice_oceandrag_coeffs.f seaice_calc_strainrates.f seaice_calc_viscosities.f seaice_calc_rhs.f seaice_calc_lhs.f
C     3. copy files seaice_*_tl.f to the corresponding seaice_*.f files, 
C        e.g. with this bash script:
C     for file in `ls seaice_*_tl.f`; do 
C       nfile=`echo $file | awk -F_ '{printf "%s_%s_%s.f", $1,$2,$3}'`; 
C       \cp -f $file $nfile
C     done
C     4. add the following code, change "call g_seaice_calc_residual" to "call seaice_calc_residual_tl", and comment out the finite 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
CMLC     Initialise
CML      DO bj=myByLo(myThid),myByHi(myThid)
CML       DO bi=myBxLo(myThid),myBxHi(myThid)
CML        DO J=1-Oly,sNy+Oly
CML         DO I=1-Olx,sNx+Olx
CML          g_duIce(I,J,bi,bj) = duice(I,J,bi,bj)
CML          g_dvIce(I,J,bi,bj) = dvice(I,J,bi,bj)
CML          g_uIceRes(I,J,bi,bj) = 0. _d 0
CML          g_vIceRes(I,J,bi,bj) = 0. _d 0
CML          uIceResP(I,J,bi,bj) = 0. _d 0
CML          vIceResP(I,J,bi,bj) = 0. _d 0
CML         ENDDO
CML        ENDDO
CML       ENDDO
CML      ENDDO
CML
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)
CML        DO J=1-Oly,sNy+Oly
CML         DO I=1-Olx,sNx+Olx
CML          duice(I,J,bi,bj)=g_uiceres(I,J,bi,bj)
CML          dvice(I,J,bi,bj)=g_viceres(I,J,bi,bj)
CML         ENDDO
CML        ENDDO
CML       ENDDO
CML      ENDDO

C     Initialise
      epsilon = SEAICE_JFNKepsilon
      reps    = 1. _d 0/epsilon

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-Oly,sNy+Oly
         DO I=1-Olx,sNx+Olx
          utp(I,J,bi,bj) = uIce(I,J,bi,bj) + epsilon * duIce(I,J,bi,bj)
          vtp(I,J,bi,bj) = vIce(I,J,bi,bj) + epsilon * dvIce(I,J,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C     Compute new residual F(u)
      CALL SEAICE_CALC_RESIDUAL(
     I     utp, vtp,
     O     uIceResP, vIceResP,
     I     newtonIter, krylovIter, myTime, myIter, myThid )

C     approximate Jacobian times vector by one-sided finite differences
C     and store in du/vIce
      DO bj = myByLo(myThid),myByHi(myThid)
       DO bi = myBxLo(myThid),myBxHi(myThid)
        DO I = 1, sNx
         DO J = 1, sNy
          duIce(I,J,bi,bj) = 
     &         (uIceResP(I,J,bi,bj)-uIceRes(I,J,bi,bj))*reps
          dvIce(I,J,bi,bj) =
     &         (vIceResP(I,J,bi,bj)-vIceRes(I,J,bi,bj))*reps
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */

      RETURN
      END
1	mlosch	1.4	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.3 2012/11/06 13:10:14 mlosch Exp $
2	mlosch	1.1	C $Name: $
3
4			#include "SEAICE_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: SEAICE_JACVEC
8			C !INTERFACE:
9			SUBROUTINE SEAICE_JACVEC(
10			I uIceLoc, vIceLoc, uIceRes, vIceRes,
11			U duIce, dvIce,
12			I newtonIter, krylovIter, myTime, myIter, myThid )
13
14			C !DESCRIPTION: \bv
15			C ==========================================================
16			C \| SUBROUTINE SEAICE_JACVEC
17			C \| o For Jacobian-free Newton-Krylov solver compute
18			C \| Jacobian times vector by finite difference approximation
19			C ==========================================================
20			C \| written by Martin Losch, Oct 2012
21			C ==========================================================
22			C \ev
23
24			C !USES:
25			IMPLICIT NONE
26
27			C === Global variables ===
28			#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31			#include "DYNVARS.h"
32			#include "GRID.h"
33			#include "SEAICE_SIZE.h"
34			#include "SEAICE_PARAMS.h"
35			#include "SEAICE.h"
36
37			C !INPUT/OUTPUT PARAMETERS:
38			C === Routine arguments ===
39			C myTime :: Simulation time
40			C myIter :: Simulation timestep number
41			C myThid :: my Thread Id. number
42			C newtonIter :: current iterate of Newton iteration
43			C krylovIter :: current iterate of Krylov iteration
44			_RL myTime
45			INTEGER myIter
46			INTEGER myThid
47			INTEGER newtonIter
48			INTEGER krylovIter
49			C u/vIceLoc :: local copies of the current ice velocity
50			_RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51			_RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
52			C u/vIceRes :: initial residual of this Newton iterate
53			_RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54			_RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55			C du/vIce :: correction of ice velocities
56			_RL duIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
57			_RL dvIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
58
59			#if ( defined (SEAICE_CGRID) && \
60			defined (SEAICE_ALLOW_JFNK) && \
61			defined (SEAICE_ALLOW_DYNAMICS) )
62			C Local variables:
63			_RL utp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
64			_RL vtp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65			C u/vIceResP :: residual computed with u/vtp
66			_RL uIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
67			_RL vIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68
69			C i,j,bi,bj :: loop indices
70			INTEGER i,j,bi,bj
71			_RL epsilon, reps
72			CEOP
73	mlosch	1.2	C Instructions for using TAF or TAMC to generate exact Jacobian times
74			C vector operations:
75			C
76			C 1. make small_f
77			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
78			C 3. staf -v1 -forward -toplevel seaice_calc_residual -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
79			C 4. insert content of taf_input_ftl.f at the end of this file
80			C 5. add the following code and comment out the finite difference code
81	mlosch	1.4	C
82			C Instruction for using TAF 2.4 and higher (or staf with default -v2
83			C starting with version 2.0):
84			C
85			C 1. make small_f
86			C 2. staf -forward -toplevel seaice_calc_residual -input uIceLoc,viceLoc -output uIceRes,vIceRes seaice_calc_residual.f seaice_oceandrag_coeffs.f seaice_calc_strainrates.f seaice_calc_viscosities.f seaice_calc_rhs.f seaice_calc_lhs.f
87			C 3. copy files seaice__tl.f to the corresponding seaice_.f files,
88			C e.g. with this bash script:
89			C for file in `ls seaice_*_tl.f`; do
90			C nfile=`echo $file \| awk -F_ '{printf "%s_%s_%s.f", $1,$2,$3}'`;
91			C \cp -f $file $nfile
92			C done
93			C 4. add the following code, change "call g_seaice_calc_residual" to "call seaice_calc_residual_tl", and comment out the finite difference code
94	mlosch	1.2	CML _RL g_duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
95			CML _RL g_dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
96			CML _RL g_uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
97			CML _RL g_vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
98			CML
99			CMLC Initialise
100			CML DO bj=myByLo(myThid),myByHi(myThid)
101			CML DO bi=myBxLo(myThid),myBxHi(myThid)
102			CML DO J=1-Oly,sNy+Oly
103			CML DO I=1-Olx,sNx+Olx
104			CML g_duIce(I,J,bi,bj) = duice(I,J,bi,bj)
105			CML g_dvIce(I,J,bi,bj) = dvice(I,J,bi,bj)
106			CML g_uIceRes(I,J,bi,bj) = 0. _d 0
107			CML g_vIceRes(I,J,bi,bj) = 0. _d 0
108			CML uIceResP(I,J,bi,bj) = 0. _d 0
109			CML vIceResP(I,J,bi,bj) = 0. _d 0
110			CML ENDDO
111			CML ENDDO
112			CML ENDDO
113			CML ENDDO
114			CML
115			CML CALL G_SEAICE_CALC_RESIDUAL( uIce, g_duice, vIce,
116			CML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
117			CML $kryloviter, mytime, myiter, mythid )
118	mlosch	1.4	CMLCML For staf -v2 replace the above with the below call
119			CMLCML CALL SEAICE_CALC_RESIDUAL_TL( uIce, g_duice, vIce,
120			CMLCML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
121			CMLCML $kryloviter, mytime, myiter, mythid )
122	mlosch	1.2	CML
123			CML DO bj=myByLo(myThid),myByHi(myThid)
124			CML DO bi=myBxLo(myThid),myBxHi(myThid)
125			CML DO J=1-Oly,sNy+Oly
126			CML DO I=1-Olx,sNx+Olx
127			CML duice(I,J,bi,bj)=g_uiceres(I,J,bi,bj)
128			CML dvice(I,J,bi,bj)=g_viceres(I,J,bi,bj)
129			CML ENDDO
130			CML ENDDO
131			CML ENDDO
132			CML ENDDO
133
134	mlosch	1.1	C Initialise
135	mlosch	1.4	epsilon = SEAICE_JFNKepsilon
136	mlosch	1.1	reps = 1. _d 0/epsilon
137
138			DO bj=myByLo(myThid),myByHi(myThid)
139			DO bi=myBxLo(myThid),myBxHi(myThid)
140			DO J=1-Oly,sNy+Oly
141			DO I=1-Olx,sNx+Olx
142			utp(I,J,bi,bj) = uIce(I,J,bi,bj) + epsilon * duIce(I,J,bi,bj)
143			vtp(I,J,bi,bj) = vIce(I,J,bi,bj) + epsilon * dvIce(I,J,bi,bj)
144			ENDDO
145			ENDDO
146			ENDDO
147			ENDDO
148
149			C Compute new residual F(u)
150			CALL SEAICE_CALC_RESIDUAL(
151			I utp, vtp,
152			O uIceResP, vIceResP,
153			I newtonIter, krylovIter, myTime, myIter, myThid )
154
155			C approximate Jacobian times vector by one-sided finite differences
156			C and store in du/vIce
157			DO bj = myByLo(myThid),myByHi(myThid)
158			DO bi = myBxLo(myThid),myBxHi(myThid)
159			DO I = 1, sNx
160			DO J = 1, sNy
161			duIce(I,J,bi,bj) =
162			& (uIceResP(I,J,bi,bj)-uIceRes(I,J,bi,bj))*reps
163			dvIce(I,J,bi,bj) =
164			& (vIceResP(I,J,bi,bj)-vIceRes(I,J,bi,bj))*reps
165			ENDDO
166			ENDDO
167			ENDDO
168			ENDDO
169
170			#endif /* SEAICE_ALLOW_DYNAMICS and SEAICE_CGRID and SEAICE_ALLOW_JFNK */
171
172			RETURN
173			END