pkg/seaice/seaice_jacvec.F

2fd913c523 Mart*0001 #include "SEAICE_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: SEAICE_JACVEC
                0005 C     !INTERFACE:
0320e25227 Mart*0006       SUBROUTINE SEAICE_JACVEC(
2fd913c523 Mart*0007      I     uIceLoc, vIceLoc, uIceRes, vIceRes,
0320e25227 Mart*0008      U     duIce, dvIce,
2fd913c523 Mart*0009      I     newtonIter, krylovIter, myTime, myIter, myThid )
                0010
                0011 C     !DESCRIPTION: \bv
                0012 C     *==========================================================*
                0013 C     | SUBROUTINE SEAICE_JACVEC
                0014 C     | o For Jacobian-free Newton-Krylov solver compute
                0015 C     |   Jacobian times vector by finite difference approximation
                0016 C     *==========================================================*
                0017 C     | written by Martin Losch, Oct 2012
                0018 C     *==========================================================*
                0019 C     \ev
                0020
                0021 C     !USES:
                0022       IMPLICIT NONE
                0023
                0024 C     === Global variables ===
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "DYNVARS.h"
                0029 #include "GRID.h"
                0030 #include "SEAICE_SIZE.h"
                0031 #include "SEAICE_PARAMS.h"
                0032 #include "SEAICE.h"
                0033
                0034 C     !INPUT/OUTPUT PARAMETERS:
                0035 C     === Routine arguments ===
                0036 C     myTime :: Simulation time
                0037 C     myIter :: Simulation timestep number
                0038 C     myThid :: my Thread Id. number
                0039 C     newtonIter :: current iterate of Newton iteration
                0040 C     krylovIter :: current iterate of Krylov iteration
                0041       _RL     myTime
                0042       INTEGER myIter
                0043       INTEGER myThid
                0044       INTEGER newtonIter
                0045       INTEGER krylovIter
                0046 C     u/vIceLoc :: local copies of the current ice velocity
                0047       _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049 C     u/vIceRes :: initial residual of this Newton iterate
                0050       _RL uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0051       _RL vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0052 C     du/vIce   :: correction of ice velocities
                0053       _RL duIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0054       _RL dvIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0055
9e7d250530 Mart*0056 #ifdef SEAICE_ALLOW_JFNK
2fd913c523 Mart*0057 C     Local variables:
                0058       _RL utp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0059       _RL vtp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0060 C     u/vIceResP :: residual computed with u/vtp
                0061       _RL uIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0062       _RL vIceResP(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0063
                0064 C     i,j,bi,bj :: loop indices
                0065       INTEGER i,j,bi,bj
                0066       _RL epsilon, reps
                0067 CEOP
cc9fc912d5 Mart*0068 C     Instructions for using TAF or TAMC to generate exact Jacobian times
70ed0fd3e4 Mart*0069 C     vector operations (if SEAICE_ALLOW_MOM_ADVECTION is defined, the
0320e25227 Mart*0070 C     file list also needs to include seaice_mom_advection.f,
                0071 C     mom_calc_hfacz.f, mom_calc_ke.f, mom_calc_relvort3.f,
70ed0fd3e4 Mart*0072 C     mom_vi_u_coriolis.f, mom_vi_u_coriolis_c4.f, mom_vi_u_grad_ke.f,
                0073 C     mom_vi_v_coriolis.f, mom_vi_v_coriolis_c4.f, mom_vi_v_grad_ke.f
                0074 C     plus flow information for diagnostics_fill.f:
                0075 CCCCCCCADJ SUBROUTINE DIAGNOSTICS_FILL INPUT  = 1,2,3,4,5,6,7,8
                0076 CCCCCCCADJ SUBROUTINE DIAGNOSTICS_FILL OUTPUT =
                0077 C     )
cc9fc912d5 Mart*0078 C
                0079 C     1. make small_f
d1e5de59de Mart*0080 C     2. cat seaice_calc_residual.f seaice_oceandrag_coeffs.f \
70ed0fd3e4 Mart*0081 C        seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f \
d1e5de59de Mart*0082 C        seaice_calc_strainrates.f seaice_calc_viscosities.f \
                0083 C        seaice_calc_rhs.f seaice_calc_lhs.f > taf_input.f
                0084 C     3. staf -v1 -forward -toplevel seaice_calc_residual \
                0085 C             -input uIceLoc,viceLoc -output uIceRes,vIceRes taf_input.f
cc9fc912d5 Mart*0086 C     4. insert content of taf_input_ftl.f at the end of this file
                0087 C     5. add the following code and comment out the finite difference code
60b28cbd4e Mart*0088 C
                0089 C     Instruction for using TAF 2.4 and higher (or staf with default -v2
                0090 C     starting with version 2.0):
                0091 C
                0092 C     1. make small_f
d1e5de59de Mart*0093 C     2. files="seaice_calc_residual.f seaice_oceandrag_coeffs.f \
70ed0fd3e4 Mart*0094 C               seaice_bottomdrag_coeffs.f seaice_calc_stressdiv.f \
d1e5de59de Mart*0095 C               seaice_calc_strainrates.f seaice_calc_viscosities.f \
                0096 C               seaice_calc_rhs.f seaice_calc_lhs.f"
                0097 C     3. staf -forward -toplevel seaice_calc_residual \
                0098 C             -input uIceLoc,viceLoc -output uIceRes,vIceRes $files
0320e25227 Mart*0099 C     4. copy files seaice_*_tl.f to the corresponding seaice_*.f files,
60b28cbd4e Mart*0100 C        e.g. with this bash script:
0320e25227 Mart*0101 C     for file in $files; do
                0102 C       nfile=`echo $file | awk -F. '{printf "%s_tl.f", $1}'`;
d1e5de59de Mart*0103 C       \cp -f $nfile $file
60b28cbd4e Mart*0104 C     done
0320e25227 Mart*0105 C     5. add the following code, change "call g_seaice_calc_residual"
                0106 C        to "call seaice_calc_residual_tl", and comment out the finite
d1e5de59de Mart*0107 C        difference code
cc9fc912d5 Mart*0108 CML      _RL g_duIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0109 CML      _RL g_dvIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0110 CML      _RL g_uIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0111 CML      _RL g_vIceRes(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0112 CML
                0113 CMLC     Initialise
                0114 CML      DO bj=myByLo(myThid),myByHi(myThid)
                0115 CML       DO bi=myBxLo(myThid),myBxHi(myThid)
                0116 CML        DO J=1-Oly,sNy+Oly
                0117 CML         DO I=1-Olx,sNx+Olx
                0118 CML          g_duIce(I,J,bi,bj) = duice(I,J,bi,bj)
                0119 CML          g_dvIce(I,J,bi,bj) = dvice(I,J,bi,bj)
                0120 CML          g_uIceRes(I,J,bi,bj) = 0. _d 0
                0121 CML          g_vIceRes(I,J,bi,bj) = 0. _d 0
                0122 CML          uIceResP(I,J,bi,bj) = 0. _d 0
                0123 CML          vIceResP(I,J,bi,bj) = 0. _d 0
                0124 CML         ENDDO
                0125 CML        ENDDO
                0126 CML       ENDDO
                0127 CML      ENDDO
                0128 CML
0320e25227 Mart*0129 CML      CALL G_SEAICE_CALC_RESIDUAL( uIce, g_duice, vIce,
                0130 CML     $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
cc9fc912d5 Mart*0131 CML     $kryloviter, mytime, myiter, mythid )
0320e25227 Mart*0132 CMLCML      For staf -v2 replace the above with the below call
                0133 CMLCML      CALL SEAICE_CALC_RESIDUAL_TL( uIce, g_duice, vIce,
                0134 CMLCML     $g_dvice, uiceresp, g_uiceres, viceresp, g_viceres, newtoniter,
60b28cbd4e Mart*0135 CMLCML     $kryloviter, mytime, myiter, mythid )
cc9fc912d5 Mart*0136 CML
                0137 CML      DO bj=myByLo(myThid),myByHi(myThid)
                0138 CML       DO bi=myBxLo(myThid),myBxHi(myThid)
                0139 CML        DO J=1-Oly,sNy+Oly
                0140 CML         DO I=1-Olx,sNx+Olx
                0141 CML          duice(I,J,bi,bj)=g_uiceres(I,J,bi,bj)
                0142 CML          dvice(I,J,bi,bj)=g_viceres(I,J,bi,bj)
                0143 CML         ENDDO
                0144 CML        ENDDO
                0145 CML       ENDDO
                0146 CML      ENDDO
                0147
2fd913c523 Mart*0148 C     Initialise
60b28cbd4e Mart*0149       epsilon = SEAICE_JFNKepsilon
2fd913c523 Mart*0150       reps    = 1. _d 0/epsilon
                0151
                0152       DO bj=myByLo(myThid),myByHi(myThid)
                0153        DO bi=myBxLo(myThid),myBxHi(myThid)
                0154         DO J=1-Oly,sNy+Oly
                0155          DO I=1-Olx,sNx+Olx
                0156           utp(I,J,bi,bj) = uIce(I,J,bi,bj) + epsilon * duIce(I,J,bi,bj)
                0157           vtp(I,J,bi,bj) = vIce(I,J,bi,bj) + epsilon * dvIce(I,J,bi,bj)
                0158          ENDDO
                0159         ENDDO
                0160        ENDDO
                0161       ENDDO
                0162
                0163 C     Compute new residual F(u)
                0164       CALL SEAICE_CALC_RESIDUAL(
                0165      I     utp, vtp,
                0166      O     uIceResP, vIceResP,
                0167      I     newtonIter, krylovIter, myTime, myIter, myThid )
                0168
                0169 C     approximate Jacobian times vector by one-sided finite differences
                0170 C     and store in du/vIce
                0171       DO bj = myByLo(myThid),myByHi(myThid)
                0172        DO bi = myBxLo(myThid),myBxHi(myThid)
                0173         DO I = 1, sNx
                0174          DO J = 1, sNy
0320e25227 Mart*0175           duIce(I,J,bi,bj) =
2fd913c523 Mart*0176      &         (uIceResP(I,J,bi,bj)-uIceRes(I,J,bi,bj))*reps
                0177           dvIce(I,J,bi,bj) =
                0178      &         (vIceResP(I,J,bi,bj)-vIceRes(I,J,bi,bj))*reps
                0179          ENDDO
                0180         ENDDO
                0181        ENDDO
                0182       ENDDO
                0183
9e7d250530 Mart*0184 #endif /* SEAICE_ALLOW_JFNK */
2fd913c523 Mart*0185
                0186       RETURN
                0187       END