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	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.3 2012/11/06 13:10:14 mlosch Exp $
2	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	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	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	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	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	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	C Initialise
135	epsilon = SEAICE_JFNKepsilon
136	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