pkg/seaice/seaice_jacvec.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.6 2013/04/23 08:33:15 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)

#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)
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 \
C        seaice_bottomdrag_coeffs.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 \
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
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_JFNK */

      RETURN
      END
1	mlosch	1.7	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.6 2013/04/23 08:33:15 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	mlosch	1.5	#ifdef SEAICE_ALLOW_JFNK
60	mlosch	1.1	C Local variables:
61			_RL utp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62			_RL vtp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
63			C u/vIceResP :: residual computed with u/vtp
64			_RL uIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
65			_RL vIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
66
67			C i,j,bi,bj :: loop indices
68			INTEGER i,j,bi,bj
69			_RL epsilon, reps
70			CEOP
71	mlosch	1.2	C Instructions for using TAF or TAMC to generate exact Jacobian times
72			C vector operations:
73			C
74			C 1. make small_f
75	mlosch	1.6	C 2. cat seaice_calc_residual.f seaice_oceandrag_coeffs.f \
76	mlosch	1.7	C seaice_bottomdrag_coeffs.f \
77	mlosch	1.6	C seaice_calc_strainrates.f seaice_calc_viscosities.f \
78			C seaice_calc_rhs.f seaice_calc_lhs.f > taf_input.f
79			C 3. staf -v1 -forward -toplevel seaice_calc_residual \
80			C -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
81	mlosch	1.2	C 4. insert content of taf_input_ftl.f at the end of this file
82			C 5. add the following code and comment out the finite difference code
83	mlosch	1.4	C
84			C Instruction for using TAF 2.4 and higher (or staf with default -v2
85			C starting with version 2.0):
86			C
87			C 1. make small_f
88	mlosch	1.6	C 2. files="seaice_calc_residual.f seaice_oceandrag_coeffs.f \
89	mlosch	1.7	C seaice_bottomdrag_coeffs.f \
90	mlosch	1.6	C seaice_calc_strainrates.f seaice_calc_viscosities.f \
91			C seaice_calc_rhs.f seaice_calc_lhs.f"
92			C 3. staf -forward -toplevel seaice_calc_residual \
93			C -input uIceLoc,viceLoc -output uIceRes,vIceRes $files
94			C 4. copy files seaice__tl.f to the corresponding seaice_.f files,
95	mlosch	1.4	C e.g. with this bash script:
96	mlosch	1.6	C for file in $files; do
97			C nfile=`echo $file \| awk -F. '{printf "%s_tl.f", $1}'`;
98			C \cp -f $nfile $file
99	mlosch	1.4	C done
100	mlosch	1.6	C 5. add the following code, change "call g_seaice_calc_residual"
101			C to "call seaice_calc_residual_tl", and comment out the finite
102			C difference code
103	mlosch	1.2	CML _RL g_duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
104			CML _RL g_dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
105			CML _RL g_uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
106			CML _RL g_vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
107			CML
108			CMLC Initialise
109			CML DO bj=myByLo(myThid),myByHi(myThid)
110			CML DO bi=myBxLo(myThid),myBxHi(myThid)
111			CML DO J=1-Oly,sNy+Oly
112			CML DO I=1-Olx,sNx+Olx
113			CML g_duIce(I,J,bi,bj) = duice(I,J,bi,bj)
114			CML g_dvIce(I,J,bi,bj) = dvice(I,J,bi,bj)
115			CML g_uIceRes(I,J,bi,bj) = 0. _d 0
116			CML g_vIceRes(I,J,bi,bj) = 0. _d 0
117			CML uIceResP(I,J,bi,bj) = 0. _d 0
118			CML vIceResP(I,J,bi,bj) = 0. _d 0
119			CML ENDDO
120			CML ENDDO
121			CML ENDDO
122			CML ENDDO
123			CML
124			CML CALL G_SEAICE_CALC_RESIDUAL( uIce, g_duice, vIce,
125			CML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
126			CML $kryloviter, mytime, myiter, mythid )
127	mlosch	1.4	CMLCML For staf -v2 replace the above with the below call
128			CMLCML CALL SEAICE_CALC_RESIDUAL_TL( uIce, g_duice, vIce,
129			CMLCML $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
130			CMLCML $kryloviter, mytime, myiter, mythid )
131	mlosch	1.2	CML
132			CML DO bj=myByLo(myThid),myByHi(myThid)
133			CML DO bi=myBxLo(myThid),myBxHi(myThid)
134			CML DO J=1-Oly,sNy+Oly
135			CML DO I=1-Olx,sNx+Olx
136			CML duice(I,J,bi,bj)=g_uiceres(I,J,bi,bj)
137			CML dvice(I,J,bi,bj)=g_viceres(I,J,bi,bj)
138			CML ENDDO
139			CML ENDDO
140			CML ENDDO
141			CML ENDDO
142
143	mlosch	1.1	C Initialise
144	mlosch	1.4	epsilon = SEAICE_JFNKepsilon
145	mlosch	1.1	reps = 1. _d 0/epsilon
146
147			DO bj=myByLo(myThid),myByHi(myThid)
148			DO bi=myBxLo(myThid),myBxHi(myThid)
149			DO J=1-Oly,sNy+Oly
150			DO I=1-Olx,sNx+Olx
151			utp(I,J,bi,bj) = uIce(I,J,bi,bj) + epsilon * duIce(I,J,bi,bj)
152			vtp(I,J,bi,bj) = vIce(I,J,bi,bj) + epsilon * dvIce(I,J,bi,bj)
153			ENDDO
154			ENDDO
155			ENDDO
156			ENDDO
157
158			C Compute new residual F(u)
159			CALL SEAICE_CALC_RESIDUAL(
160			I utp, vtp,
161			O uIceResP, vIceResP,
162			I newtonIter, krylovIter, myTime, myIter, myThid )
163
164			C approximate Jacobian times vector by one-sided finite differences
165			C and store in du/vIce
166			DO bj = myByLo(myThid),myByHi(myThid)
167			DO bi = myBxLo(myThid),myBxHi(myThid)
168			DO I = 1, sNx
169			DO J = 1, sNy
170			duIce(I,J,bi,bj) =
171			& (uIceResP(I,J,bi,bj)-uIceRes(I,J,bi,bj))*reps
172			dvIce(I,J,bi,bj) =
173			& (vIceResP(I,J,bi,bj)-vIceRes(I,J,bi,bj))*reps
174			ENDDO
175			ENDDO
176			ENDDO
177			ENDDO
178
179	mlosch	1.5	#endif /* SEAICE_ALLOW_JFNK */
180	mlosch	1.1
181			RETURN
182			END