C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_jacvec.F,v 1.8 2017/05/23 16:24:46 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 (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 \ 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 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