C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/dgoldberg/streamice/adstreamice_cg_solve.F,v 1.1 2012/07/19 18:52:31 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 iters, I myThid ) C /============================================================\ C | SUBROUTINE | C | o | C |============================================================| C | | C \============================================================/ 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 iters 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 ! iters = streamice_max_cg_iter #ifdef ALLOW_STREAMICE PRINT *, "CALLING MANUAL CG ADJOINT" 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 CALL STREAMICE_CG_SOLVE( & U_state, & V_state, & Bu_state, & Bv_state, & A_uu, & A_uv, & A_vu, & A_vv, & tolerance, & tmpiter, & myThid ) tmpiter = 0 CALL STREAMICE_CG_SOLVE( & cg_Uin, & cg_Vin, & cg_Bu, & cg_Bv, & A_uu, & A_uv, & A_vu, & A_vv, & ad_tolerance, & tmpiter, & myThid ) ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1,sNy ! DO i=1,sNx ! Zu_SI (i,j,bi,bj) = 0. _d 0 ! Zv_SI (i,j,bi,bj) = 0. _d 0 ! Ru_SI (i,j,bi,bj) = 0. _d 0 ! Rv_SI (i,j,bi,bj) = 0. _d 0 ! Au_SI (i,j,bi,bj) = 0. _d 0 ! Av_SI (i,j,bi,bj) = 0. _d 0 ! Du_SI (i,j,bi,bj) = 0. _d 0 ! Dv_SI (i,j,bi,bj) = 0. _d 0 ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! C FIND INITIAL RESIDUAL, and initialize r ! ! ! #ifdef STREAMICE_CONSTRUCT_MATRIX ! ! ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=js,je ! DO i=is,ie ! DO colx=-1,1 ! DO coly=-1,1 ! Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) + ! & A_uu(i,j,bi,bj,colx,coly)* ! & cg_Uin(i+colx,j+coly,bi,bj)+ ! & A_uv(i,j,bi,bj,colx,coly)* ! & cg_Vin(i+colx,j+coly,bi,bj) ! Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) + ! & A_vu(i,j,bi,bj,colx,coly)* ! & cg_Uin(i+colx,j+coly,bi,bj)+ ! & A_vv(i,j,bi,bj,colx,coly)* ! & cg_Vin(i+colx,j+coly,bi,bj) ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! ! ! print *, "GOT HERE 2" ! ! ! #else ! ! ! ! CALL STREAMICE_CG_ACTION( myThid, ! ! O Au_SI, ! ! O Av_SI, ! ! I cg_Uin, ! ! I cg_Vin, ! ! I 0, sNx+1, 0, sNy+1 ) ! ! ! ! #endif ! ! _EXCH_XY_RL( Au_SI, myThid ) ! _EXCH_XY_RL( Av_SI, myThid ) ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1-OLy,sNy+OLy ! DO i=1-OLx,sNx+OLx ! Ru_SI(i,j,bi,bj)=cg_Bu(i,j,bi,bj)- ! & Au_SI(i,j,bi,bj) ! Rv_SI(i,j,bi,bj)=cg_Bv(i,j,bi,bj)- ! & Av_SI(i,j,bi,bj) ! ENDDO ! ENDDO ! dot_p1_tile(bi,bj) = 0. _d 0 ! dot_p2_tile(bi,bj) = 0. _d 0 ! ENDDO ! ENDDO ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1,sNy ! DO i=1,sNx ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2 ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2 ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) ! resid_0 = sqrt(dot_p1) ! ! C CCCCCCCCCCCCCCCCCCCC ! ! 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 (STREAMICE_umask(i,j,bi,bj).eq.1.0) ! & Zu_SI(i,j,bi,bj)=Ru_SI(i,j,bi,bj) / DIAGu_SI(i,j,bi,bj) ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) ! & Zv_SI(i,j,bi,bj)=Rv_SI(i,j,bi,bj) / DIAGv_SI(i,j,bi,bj) ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! cg_halo = min(OLx-1,OLy-1) ! conv_flag = 0 ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1-OLy,sNy+OLy ! DO i=1-OLx,sNx+OLx ! Du_SI(i,j,bi,bj)=Zu_SI(i,j,bi,bj) ! Dv_SI(i,j,bi,bj)=Zv_SI(i,j,bi,bj) ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! resid = resid_0 ! iters = 0 ! ! c !!!!!!!!!!!!!!!!!! ! c !! !! ! c !! MAIN CG LOOP !! ! c !! !! ! c !!!!!!!!!!!!!!!!!! ! ! ! ! ! c ! initially, b-grid data is valid up to 3 halo nodes out -- right? (check for MITgcm!!) ! ! ! ! ! IF(STREAMICE_construct_matrix) CALL STREAMICE_CG_MAKE_A(myThid) ! ! ! do iter = 1, streamice_max_cg_iter ! if (resid .gt. tolerance*resid_0) then ! ! c to avoid using "exit" ! iters = iters + 1 ! ! is = 1 - cg_halo ! ie = sNx + cg_halo ! js = 1 - cg_halo ! je = sNy + cg_halo ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1-OLy,sNy+OLy ! DO i=1-OLx,sNx+OLx ! Au_SI(i,j,bi,bj) = 0. _d 0 ! Av_SI(i,j,bi,bj) = 0. _d 0 ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! ! IF (STREAMICE_construct_matrix) THEN ! ! ! #ifdef STREAMICE_CONSTRUCT_MATRIX ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=js,je ! DO i=is,ie ! DO colx=-1,1 ! DO coly=-1,1 ! Au_SI(i,j,bi,bj) = Au_SI(i,j,bi,bj) + ! & A_uu(i,j,bi,bj,colx,coly)* ! & Du_SI(i+colx,j+coly,bi,bj)+ ! & A_uv(i,j,bi,bj,colx,coly)* ! & Dv_SI(i+colx,j+coly,bi,bj) ! Av_SI(i,j,bi,bj) = Av_SI(i,j,bi,bj) + ! & A_vu(i,j,bi,bj,colx,coly)* ! & Du_SI(i+colx,j+coly,bi,bj)+ ! & A_vv(i,j,bi,bj,colx,coly)* ! & Dv_SI(i+colx,j+coly,bi,bj) ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! ! else ! ! #else ! ! ! ! CALL STREAMICE_CG_ACTION( myThid, ! ! O Au_SI, ! ! O Av_SI, ! ! I Du_SI, ! ! I Dv_SI, ! ! I is,ie,js,je) ! ! ! ! ! ENDIF ! ! ! ! #endif ! ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! dot_p1_tile(bi,bj) = 0. _d 0 ! dot_p2_tile(bi,bj) = 0. _d 0 ! ENDDO ! ENDDO ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1,sNy ! DO i=1,sNx ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)* ! & Ru_SI(i,j,bi,bj) ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Du_SI(i,j,bi,bj)* ! & Au_SI(i,j,bi,bj) ! ENDIF ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)* ! & Rv_SI(i,j,bi,bj) ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Dv_SI(i,j,bi,bj)* ! & Av_SI(i,j,bi,bj) ! ENDIF ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) ! CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid ) ! alpha_k = dot_p1/dot_p2 ! ! 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 (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN ! cg_Uin(i,j,bi,bj)=cg_Uin(i,j,bi,bj)+ ! & alpha_k*Du_SI(i,j,bi,bj) ! Ru_old_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) ! Zu_old_SI(i,j,bi,bj) = Zu_SI(i,j,bi,bj) ! Ru_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj)- ! & alpha_k*Au_SI(i,j,bi,bj) ! Zu_SI(i,j,bi,bj) = Ru_SI(i,j,bi,bj) / ! & DIAGu_SI(i,j,bi,bj) ! ENDIF ! ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN ! cg_Vin(i,j,bi,bj)=cg_Vin(i,j,bi,bj)+ ! & alpha_k*Dv_SI(i,j,bi,bj) ! Rv_old_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) ! Zv_old_SI(i,j,bi,bj) = Zv_SI(i,j,bi,bj) ! Rv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj)- ! & alpha_k*Av_SI(i,j,bi,bj) ! Zv_SI(i,j,bi,bj) = Rv_SI(i,j,bi,bj) / ! & DIAGv_SI(i,j,bi,bj) ! ! ENDIF ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! dot_p1_tile(bi,bj) = 0. _d 0 ! dot_p2_tile(bi,bj) = 0. _d 0 ! ENDDO ! ENDDO ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1,sNy ! DO i=1,sNx ! ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) THEN ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zu_SI(i,j,bi,bj)* ! & Ru_SI(i,j,bi,bj) ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zu_old_SI(i,j,bi,bj)* ! & Ru_old_SI(i,j,bi,bj) ! ENDIF ! ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) THEN ! dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Zv_SI(i,j,bi,bj)* ! & Rv_SI(i,j,bi,bj) ! dot_p2_tile(bi,bj)=dot_p2_tile(bi,bj)+Zv_old_SI(i,j,bi,bj)* ! & Rv_old_SI(i,j,bi,bj) ! ENDIF ! ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) ! CALL GLOBAL_SUM_TILE_RL( dot_p2_tile, dot_p2, myThid ) ! ! beta_k = dot_p1/dot_p2 ! ! 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 (STREAMICE_umask(i,j,bi,bj).eq.1.0) ! & Du_SI(i,j,bi,bj)=beta_k*Du_SI(i,j,bi,bj)+ ! & Zu_SI(i,j,bi,bj) ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) ! & Dv_SI(i,j,bi,bj)=beta_k*Dv_SI(i,j,bi,bj)+ ! & Zv_SI(i,j,bi,bj) ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! dot_p1_tile(bi,bj) = 0. _d 0 ! ENDDO ! ENDDO ! ! DO bj = myByLo(myThid), myByHi(myThid) ! DO bi = myBxLo(myThid), myBxHi(myThid) ! DO j=1,sNy ! DO i=1,sNx ! IF (STREAMICE_umask(i,j,bi,bj).eq.1.0) ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Ru_SI(i,j,bi,bj)**2 ! IF (STREAMICE_vmask(i,j,bi,bj).eq.1.0) ! & dot_p1_tile(bi,bj)=dot_p1_tile(bi,bj)+Rv_SI(i,j,bi,bj)**2 ! ENDDO ! ENDDO ! ENDDO ! ENDDO ! ! CALL GLOBAL_SUM_TILE_RL( dot_p1_tile, dot_p1, myThid ) ! resid = sqrt(dot_p1) ! ! ! IF (iter .eq. 1) then ! ! print *, alpha_k, beta_k, resid ! ! ENDIF ! ! cg_halo = cg_halo - 1 ! ! if (cg_halo .eq. 0) then ! cg_halo = min(OLx-1,OLy-1) ! _EXCH_XY_RL( Du_SI, myThid ) ! _EXCH_XY_RL( Dv_SI, myThid ) ! _EXCH_XY_RL( Ru_SI, myThid ) ! _EXCH_XY_RL( Rv_SI, myThid ) ! _EXCH_XY_RL( cg_Uin, myThid ) ! _EXCH_XY_RL( cg_Vin, myThid ) ! endif ! ! endif ! enddo ! end of CG loop c to avoid using "exit" c if iters has reached max_iters there is no convergence ! IF (iters .lt. streamice_max_cg_iter) THEN ! conv_flag = 1 ! ENDIF DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO j=1-Oly,sNy+Oly DO i=1-Olx,sNx+Olx 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) ! print *,"ADA",cg_Uin(i,j,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 ! if (i.eq.1.and.j.eq.1) then ! print *, "adA_uu", adA_uu(i,j,bi,bj,-1,-1), ! & adA_uu(i,j,bi,bj,-1,0), ! & adA_uu(i,j,bi,bj,-1,1), ! & cg_Uin(i,j,bi,bj) ! endif enddo enddo enddo enddo ! print *,"MATRIX 1" ! do i=1,sNx ! do j=1,sNy ! print *, ada_uu(i,j,1,1,-1,-1), ada_uu(i,j,1,1,0,-1), ! & ada_uu(i,j,1,1,1,-1), ada_uu(i,j,1,1,-1,0), ! & ada_uu(i,j,1,1,0,0), ada_uu(i,j,1,1,1,0), ! & ada_uu(i,j,1,1,-1,1), ada_uu(i,j,1,1,0,1),ada_uu(i,j,1,1,1,1) ! enddo ! enddo ! ! print *,"MATRIX 2" ! do i=1,sNx ! do j=1,sNy ! print *, ada_uv(i,j,1,1,-1,-1), ada_uv(i,j,1,1,0,-1), ! & ada_uv(i,j,1,1,1,-1), ada_uv(i,j,1,1,-1,0), ! & ada_uv(i,j,1,1,0,0), ada_uv(i,j,1,1,1,0), ! & ada_uv(i,j,1,1,-1,1), ada_uv(i,j,1,1,0,1),ada_uv(i,j,1,1,1,1) ! enddo ! enddo ! ! print *,"MATRIX 3" ! do i=1,sNx ! do j=1,sNy ! print *, ada_vu(i,j,1,1,-1,-1), ada_vu(i,j,1,1,0,-1), ! & ada_vu(i,j,1,1,1,-1), ada_vu(i,j,1,1,-1,0), ! & ada_vu(i,j,1,1,0,0), ada_vu(i,j,1,1,1,0), ! & ada_vu(i,j,1,1,-1,1), ada_vu(i,j,1,1,0,1),ada_vu(i,j,1,1,1,1) ! enddo ! enddo ! ! print *,"MATRIX 4" ! do i=1,sNx ! do j=1,sNy ! print *, ada_vv(i,j,1,1,-1,-1), ada_vv(i,j,1,1,0,-1), ! & ada_vv(i,j,1,1,1,-1), ada_vv(i,j,1,1,-1,0), ! & ada_vv(i,j,1,1,0,0), ada_vv(i,j,1,1,1,0), ! & ada_vv(i,j,1,1,-1,1), ada_vv(i,j,1,1,0,1),ada_vv(i,j,1,1,1,1) ! 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 PRINT *, "DONE WITH MANUAL CG ADJOINT:",tmpiter,"ITERS" #endif RETURN END