C $Header: /u/gcmpack/MITgcm/pkg/streamice/adstreamice_cg_solve.F,v 1.7 2014/12/13 14:13:11 dgoldberg Exp $ C $Name: $ #include "STREAMICE_OPTIONS.h" C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP SUBROUTINE ADSTREAMICE_CG_SOLVE( U U_state, ! velocities - need to be recalc ed I cg_Bu, ! adjoint of vel (input) U V_state, ! velocities - need to be recalc ed I cg_Bv, ! adjoint of vel (input) I Bu_state, ! to recalc velocities U cg_Uin, ! adjoint of RHS (output) I Bv_state, ! to recalc velocities U cg_Vin, ! adjoint of RHS (output) I A_uu, ! section of matrix that multiplies u and projects on u U adA_uu, ! adjoint of matrix coeffs (output) I A_uv, ! section of matrix that multiplies v and projects on u U adA_uv, ! adjoint of matrix coeffs (output) I A_vu, ! section of matrix that multiplies u and projects on v U adA_vu, ! adjoint of matrix coeffs (output) I A_vv, ! section of matrix that multiplies v and projects on v U adA_vv, ! adjoint of matrix coeffs (output) I tolerance, I maxiters, I myThid ) C *============================================================* C | SUBROUTINE | C | o | C *============================================================* C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "STREAMICE.h" #include "STREAMICE_CG.h" C !INPUT/OUTPUT ARGUMENTS C cg_Uin, cg_Vin - input and output velocities C cg_Bu, cg_Bv - driving stress _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL U_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & adA_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & adA_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1), & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1) _RL tolerance INTEGER maxiters INTEGER myThid C !LOCAL VARIABLES INTEGER i, j, bi, bj, cg_halo, conv_flag, tmpiter INTEGER iter, is, js, ie, je, colx, coly, k _RL Utemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL Vtemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL UtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL VtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL dot_p1, dot_p2, alpha_k, beta_k, resid, resid_0 _RL dot_p1_tile (nSx,nSy) _RL dot_p2_tile (nSx,nSy) _RL ad_tolerance CHARACTER*(MAX_LEN_MBUF) msgBuf CEOP ! iters = streamice_max_cg_iter #ifdef ALLOW_STREAMICE WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , 1) ! print *, "GOT HERE mythid=", mythid, tolerance conv_flag = 0 ad_tolerance = 1.e-14 DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx Utemp (i,j,bi,bj) = & cg_Uin (i,j,bi,bj) Vtemp (i,j,bi,bj) = & cg_Vin (i,j,bi,bj) UtempSt (i,j,bi,bj) = & U_state (i,j,bi,bj) VtempSt (i,j,bi,bj) = & V_state (i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO ! print *, "GOT HERE 2 mythid=", mythid, tolerance #ifndef ALLOW_OPENAD CALL STREAMICE_CG_SOLVE( & U_state, & V_state, & Bu_state, & Bv_state, & A_uu, & A_uv, & A_vu, & A_vv, & tolerance, & tmpiter, & maxiters, & myThid ) #endif ! print *, "GOT HERE 3 mythid=", mythid, tolerance tmpiter = 0 _EXCH_XY_RL( cg_Bu, myThid ) _EXCH_XY_RL( cg_Bv, myThid ) CALL STREAMICE_CG_SOLVE( & cg_Uin, & cg_Vin, & cg_Bu, & cg_Bv, & A_uu, & A_uv, & A_vu, & A_vv, & ad_tolerance, & tmpiter, & maxiters, & myThid ) ! print *, "GOT HERE 4 mythid=", mythid, tolerance _EXCH_XY_RL( cg_Uin, myThid ) _EXCH_XY_RL( cg_Vin, myThid ) _EXCH_XY_RL( U_state, myThid ) _EXCH_XY_RL( V_state, myThid ) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1,sNy DO i=1,sNx DO colx=-1,1 DO coly=-1,1 if (STREAMICE_umask(i,j,bi,bj).eq.1) then if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then adA_uu(i,j,bi,bj,colx,coly) = & adA_uu(i,j,bi,bj,colx,coly) - & cg_Uin(i,j,bi,bj) * & U_state(i+colx,j+coly,bi,bj) endif if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then adA_uv(i,j,bi,bj,colx,coly) = & adA_uv(i,j,bi,bj,colx,coly) - & cg_Uin(i,j,bi,bj) * & V_state(i+colx,j+coly,bi,bj) endif endif if (STREAMICE_vmask(i,j,bi,bj).eq.1) then if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then adA_vu(i,j,bi,bj,colx,coly) = & adA_vu(i,j,bi,bj,colx,coly) - & cg_Vin(i,j,bi,bj) * & U_state(i+colx,j+coly,bi,bj) endif if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then adA_vv(i,j,bi,bj,colx,coly) = & adA_vv(i,j,bi,bj,colx,coly) - & cg_Vin(i,j,bi,bj) * & V_state(i+colx,j+coly,bi,bj) endif endif enddo enddo enddo enddo enddo enddo DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx if (i.lt.1.or.i.gt.sNx.or. & j.lt.1.or.j.gt.sNy) then cg_Uin (i,j,bi,bj) = 0.0 cg_Vin (i,j,bi,bj) = 0.0 DO colx=-1,1 DO coly=-1,1 ada_uu(i,j,bi,bj,colx,coly)=0.0 ada_uv(i,j,bi,bj,colx,coly)=0.0 ada_vu(i,j,bi,bj,colx,coly)=0.0 ada_vv(i,j,bi,bj,colx,coly)=0.0 enddo enddo endif cg_Uin (i,j,bi,bj) = & cg_Uin (i,j,bi,bj) + & Utemp (i,j,bi,bj) cg_Vin (i,j,bi,bj) = & cg_Vin (i,j,bi,bj) + & Vtemp (i,j,bi,bj) cg_bu (i,j,bi,bj) = 0. cg_bv (i,j,bi,bj) = 0. U_state (i,j,bi,bj) = & UtempSt (i,j,bi,bj) V_state (i,j,bi,bj) = & VtempSt (i,j,bi,bj) ENDDO ENDDO ENDDO ENDDO WRITE(msgBuf,'(A,I5,A)') 'DONE WITH MANUAL CG ADJOINT:',tmpiter, & 'ITERS' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , 1) #endif RETURN END