pkg/streamice/adstreamice_cg_solve.F

5ca83cd8f7 Dani*0001 #include "STREAMICE_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004
                0005 CBOP
f991e74a3b Jean*0006       SUBROUTINE ADSTREAMICE_CG_SOLVE(
c8fca1659b Jean*0007      U                               U_state,    ! velocities - need to be recalc ed
5ca83cd8f7 Dani*0008      I                               cg_Bu,      ! adjoint of vel (input)
c8fca1659b Jean*0009      U                               V_state,    ! velocities - need to be recalc ed
5ca83cd8f7 Dani*0010      I                               cg_Bv,      ! adjoint of vel (input)
                0011      I                               Bu_state,   ! to recalc velocities
                0012      U                               cg_Uin,     ! adjoint of RHS (output)
                0013      I                               Bv_state,   ! to recalc velocities
                0014      U                               cg_Vin,     ! adjoint of RHS (output)
                0015      I                               A_uu,       ! section of matrix that multiplies u and projects on u
                0016      U                               adA_uu,     ! adjoint of matrix coeffs (output)
                0017      I                               A_uv,       ! section of matrix that multiplies v and projects on u
                0018      U                               adA_uv,     ! adjoint of matrix coeffs (output)
                0019      I                               A_vu,       ! section of matrix that multiplies u and projects on v
                0020      U                               adA_vu,     ! adjoint of matrix coeffs (output)
                0021      I                               A_vv,       ! section of matrix that multiplies v and projects on v
                0022      U                               adA_vv,     ! adjoint of matrix coeffs (output)
f991e74a3b Jean*0023      I                               tolerance,
d2cdb9260d Dani*0024      I                               maxiters,
feb7fa5d1e dngo*0025      I                               myThid)
                0026 C     *============================================================*
                0027 C     | SUBROUTINE ADSTREAMICE_CG_SOLVE                            |
                0028 C     | o wrapper for ADSTREAMICE_CG_SOLVE_RECALC                  |
                0029 C     *============================================================*
                0030
                0031 C     !USES:
                0032       IMPLICIT NONE
                0033
                0034 C     === Global variables ===
                0035 #include "SIZE.h"
                0036 #include "EEPARAMS.h"
                0037 #include "PARAMS.h"
                0038
                0039 C     !INPUT/OUTPUT ARGUMENTS
                0040 C     cg_Uin, cg_Vin - input and output velocities
                0041 C     cg_Bu, cg_Bv - driving stress
                0042       _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0043       _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0044       _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0045       _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0046       _RL U_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0047       _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0049       _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0050       _RL
                0051      & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0052      & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0053      & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0054      & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0055      & adA_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0056      & adA_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0057      & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0058      & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
                0059       _RL tolerance
                0060       INTEGER maxiters
                0061       INTEGER myThid
                0062
                0063 #ifdef ALLOW_STREAMICE
                0064       INTEGER recalc_vel
                0065
                0066       recalc_vel = 0
                0067       CALL ADSTREAMICE_CG_SOLVE_RECALC(
                0068      U                               U_state,    ! velocities - need to be recalc ed
                0069      I                               cg_Bu,      ! adjoint of vel (input)
                0070      U                               V_state,    ! velocities - need to be recalc ed
                0071      I                               cg_Bv,      ! adjoint of vel (input)
                0072      I                               Bu_state,   ! to recalc velocities
                0073      U                               cg_Uin,     ! adjoint of RHS (output)
                0074      I                               Bv_state,   ! to recalc velocities
                0075      U                               cg_Vin,     ! adjoint of RHS (output)
                0076      I                               A_uu,       ! section of matrix that multiplies u and projects on u
                0077      U                               adA_uu,     ! adjoint of matrix coeffs (output)
                0078      I                               A_uv,       ! section of matrix that multiplies v and projects on u
                0079      U                               adA_uv,     ! adjoint of matrix coeffs (output)
                0080      I                               A_vu,       ! section of matrix that multiplies u and projects on v
                0081      U                               adA_vu,     ! adjoint of matrix coeffs (output)
                0082      I                               A_vv,       ! section of matrix that multiplies v and projects on v
                0083      U                               adA_vv,     ! adjoint of matrix coeffs (output)
                0084      I                               tolerance,
                0085      I                               maxiters,
                0086      I                               recalc_vel,
                0087      I                               myThid)
                0088
                0089 c       iters = streamice_max_cg_iter
                0090
                0091 #endif
                0092       RETURN
                0093       END
                0094 CEOP
                0095
                0096 CBOP
                0097       SUBROUTINE ADSTREAMICE_CG_SOLVE_RECALC(
                0098      U                               U_state,    ! velocities - need to be recalc ed
                0099      I                               cg_Bu,      ! adjoint of vel (input)
                0100      U                               V_state,    ! velocities - need to be recalc ed
                0101      I                               cg_Bv,      ! adjoint of vel (input)
                0102      I                               Bu_state,   ! to recalc velocities
                0103      U                               cg_Uin,     ! adjoint of RHS (output)
                0104      I                               Bv_state,   ! to recalc velocities
                0105      U                               cg_Vin,     ! adjoint of RHS (output)
                0106      I                               A_uu,       ! section of matrix that multiplies u and projects on u
                0107      U                               adA_uu,     ! adjoint of matrix coeffs (output)
                0108      I                               A_uv,       ! section of matrix that multiplies v and projects on u
                0109      U                               adA_uv,     ! adjoint of matrix coeffs (output)
                0110      I                               A_vu,       ! section of matrix that multiplies u and projects on v
                0111      U                               adA_vu,     ! adjoint of matrix coeffs (output)
                0112      I                               A_vv,       ! section of matrix that multiplies v and projects on v
                0113      U                               adA_vv,     ! adjoint of matrix coeffs (output)
                0114      I                               tolerance,
                0115      I                               maxiters,
                0116      I                               recalc_vel,
5ca83cd8f7 Dani*0117      I                               myThid )
f991e74a3b Jean*0118 C     *============================================================*
feb7fa5d1e dngo*0119 C     | SUBROUTINE ADSTREAMICE_CG_SOLVE_RECALC                     |
                0120 C     | o called by ADSTREAMICE_CG_SOLVE or directly               |
                0121 C     | o can enable recalc of velocity if needed                  |
f991e74a3b Jean*0122 C     *============================================================*
                0123
                0124 C     !USES:
5ca83cd8f7 Dani*0125       IMPLICIT NONE
                0126
                0127 C     === Global variables ===
                0128 #include "SIZE.h"
                0129 #include "EEPARAMS.h"
                0130 #include "PARAMS.h"
                0131 #include "STREAMICE.h"
                0132 #include "STREAMICE_CG.h"
                0133
                0134 C     !INPUT/OUTPUT ARGUMENTS
                0135 C     cg_Uin, cg_Vin - input and output velocities
                0136 C     cg_Bu, cg_Bv - driving stress
                0137       _RL cg_Uin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0138       _RL cg_Vin (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0139       _RL cg_Bu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0140       _RL cg_Bv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0141       _RL U_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0142       _RL V_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0143       _RL Bu_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0144       _RL Bv_state (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f991e74a3b Jean*0145       _RL
5ca83cd8f7 Dani*0146      & A_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0147      & A_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0148      & A_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0149      & A_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0150      & adA_uu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0151      & adA_vu (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0152      & adA_uv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1),
                0153      & adA_vv (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy,-1:1,-1:1)
                0154       _RL tolerance
d2cdb9260d Dani*0155       INTEGER maxiters
feb7fa5d1e dngo*0156       INTEGER recalc_vel ! 1 = recalculate forward velocity
                0157                          !     before solving adjoint
                0158                          ! 0 = velocity already correct
5ca83cd8f7 Dani*0159       INTEGER myThid
                0160
f991e74a3b Jean*0161 C     !LOCAL VARIABLES
96b006450c dngo*0162       INTEGER i, j, bi, bj, conv_flag, tmpiter
                0163       INTEGER colx, coly
5ca83cd8f7 Dani*0164       _RL Utemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0165       _RL Vtemp (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0166       _RL UtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0167       _RL VtempSt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0168       _RL ad_tolerance
                0169       CHARACTER*(MAX_LEN_MBUF) msgBuf
f991e74a3b Jean*0170 CEOP
5ca83cd8f7 Dani*0171
96b006450c dngo*0172 c       iters = streamice_max_cg_iter
5ca83cd8f7 Dani*0173
                0174 #ifdef ALLOW_STREAMICE
                0175
                0176       WRITE(msgBuf,'(A)') 'CALLING MANUAL CG ADJOINT'
                0177        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0178      &                     SQUEEZE_RIGHT , 1)
                0179
                0180       conv_flag = 0
                0181       ad_tolerance = 1.e-14
                0182
                0183       DO bj = myByLo(myThid), myByHi(myThid)
                0184        DO bi = myBxLo(myThid), myBxHi(myThid)
f991e74a3b Jean*0185         DO j=1-OLy,sNy+OLy
                0186          DO i=1-OLx,sNx+OLx
5ca83cd8f7 Dani*0187           Utemp (i,j,bi,bj) =
                0188      &     cg_Uin (i,j,bi,bj)
                0189           Vtemp (i,j,bi,bj) =
                0190      &     cg_Vin (i,j,bi,bj)
                0191           UtempSt (i,j,bi,bj) =
                0192      &     U_state (i,j,bi,bj)
                0193           VtempSt (i,j,bi,bj) =
                0194      &     V_state (i,j,bi,bj)
                0195          ENDDO
                0196         ENDDO
                0197        ENDDO
                0198       ENDDO
                0199
feb7fa5d1e dngo*0200 !#if !defined(ALLOW_OPENAD) & !defined(ALLOW_TAPENADE)
                0201       if (recalc_vel .eq. 1) then
                0202        CALL STREAMICE_CG_SOLVE(
5ca83cd8f7 Dani*0203      &  U_state,
                0204      &  V_state,
                0205      &  Bu_state,
                0206      &  Bv_state,
                0207      &  A_uu,
                0208      &  A_uv,
                0209      &  A_vu,
f991e74a3b Jean*0210      &  A_vv,
                0211      &  tolerance,
5ca83cd8f7 Dani*0212      &  tmpiter,
d2cdb9260d Dani*0213      &  maxiters,
5ca83cd8f7 Dani*0214      &  myThid )
feb7fa5d1e dngo*0215       endif
                0216 !#endif
d2cdb9260d Dani*0217
5ca83cd8f7 Dani*0218       tmpiter = 0
                0219
                0220       _EXCH_XY_RL( cg_Bu, myThid )
                0221       _EXCH_XY_RL( cg_Bv, myThid )
                0222
f991e74a3b Jean*0223       CALL STREAMICE_CG_SOLVE(
5ca83cd8f7 Dani*0224      &  cg_Uin,
                0225      &  cg_Vin,
                0226      &  cg_Bu,
                0227      &  cg_Bv,
                0228      &  A_uu,
                0229      &  A_uv,
                0230      &  A_vu,
f991e74a3b Jean*0231      &  A_vv,
                0232      &  ad_tolerance,
5ca83cd8f7 Dani*0233      &  tmpiter,
d2cdb9260d Dani*0234      &  maxiters,
5ca83cd8f7 Dani*0235      &  myThid )
                0236
                0237       _EXCH_XY_RL( cg_Uin, myThid )
f991e74a3b Jean*0238       _EXCH_XY_RL( cg_Vin, myThid )
                0239       _EXCH_XY_RL( U_state, myThid )
                0240       _EXCH_XY_RL( V_state, myThid )
5ca83cd8f7 Dani*0241
                0242       DO bj = myByLo(myThid), myByHi(myThid)
                0243        DO bi = myBxLo(myThid), myBxHi(myThid)
5721e29731 Dani*0244         DO j=1,sNy
                0245          DO i=1,sNx
5ca83cd8f7 Dani*0246           DO colx=-1,1
                0247            DO coly=-1,1
                0248
                0249             if (STREAMICE_umask(i,j,bi,bj).eq.1) then
                0250              if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
f991e74a3b Jean*0251                 adA_uu(i,j,bi,bj,colx,coly) =
                0252      &           adA_uu(i,j,bi,bj,colx,coly) -
5ca83cd8f7 Dani*0253      &           cg_Uin(i,j,bi,bj) *
                0254      &           U_state(i+colx,j+coly,bi,bj)
                0255
                0256              endif
                0257              if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
                0258                 adA_uv(i,j,bi,bj,colx,coly) =
f991e74a3b Jean*0259      &           adA_uv(i,j,bi,bj,colx,coly) -
5ca83cd8f7 Dani*0260      &           cg_Uin(i,j,bi,bj) *
                0261      &           V_state(i+colx,j+coly,bi,bj)
                0262              endif
                0263             endif
                0264
                0265             if (STREAMICE_vmask(i,j,bi,bj).eq.1) then
                0266              if (STREAMICE_umask(i+colx,j+coly,bi,bj).eq.1) then
                0267                 adA_vu(i,j,bi,bj,colx,coly) =
f991e74a3b Jean*0268      &           adA_vu(i,j,bi,bj,colx,coly) -
5ca83cd8f7 Dani*0269      &           cg_Vin(i,j,bi,bj) *
                0270      &           U_state(i+colx,j+coly,bi,bj)
                0271              endif
                0272              if (STREAMICE_vmask(i+colx,j+coly,bi,bj).eq.1) then
                0273                 adA_vv(i,j,bi,bj,colx,coly) =
f991e74a3b Jean*0274      &           adA_vv(i,j,bi,bj,colx,coly) -
5ca83cd8f7 Dani*0275      &           cg_Vin(i,j,bi,bj) *
                0276      &           V_state(i+colx,j+coly,bi,bj)
                0277              endif
                0278             endif
                0279
                0280            enddo
                0281           enddo
                0282          enddo
                0283         enddo
                0284        enddo
                0285       enddo
                0286
                0287       DO bj = myByLo(myThid), myByHi(myThid)
                0288        DO bi = myBxLo(myThid), myBxHi(myThid)
f991e74a3b Jean*0289         DO j=1-OLy,sNy+OLy
                0290          DO i=1-OLx,sNx+OLx
5ca83cd8f7 Dani*0291           if (i.lt.1.or.i.gt.sNx.or.
                0292      &        j.lt.1.or.j.gt.sNy) then
                0293            cg_Uin (i,j,bi,bj) = 0.0
                0294            cg_Vin (i,j,bi,bj) = 0.0
f991e74a3b Jean*0295
5ca83cd8f7 Dani*0296            DO colx=-1,1
f991e74a3b Jean*0297             DO coly=-1,1
5ca83cd8f7 Dani*0298              ada_uu(i,j,bi,bj,colx,coly)=0.0
                0299              ada_uv(i,j,bi,bj,colx,coly)=0.0
                0300              ada_vu(i,j,bi,bj,colx,coly)=0.0
                0301              ada_vv(i,j,bi,bj,colx,coly)=0.0
                0302             enddo
                0303            enddo
                0304
                0305           endif
                0306           cg_Uin (i,j,bi,bj) =
f991e74a3b Jean*0307      &     cg_Uin (i,j,bi,bj) +
5ca83cd8f7 Dani*0308      &     Utemp (i,j,bi,bj)
                0309           cg_Vin (i,j,bi,bj) =
f991e74a3b Jean*0310      &     cg_Vin (i,j,bi,bj) +
5ca83cd8f7 Dani*0311      &     Vtemp (i,j,bi,bj)
                0312           cg_bu (i,j,bi,bj) = 0.
                0313           cg_bv (i,j,bi,bj) = 0.
                0314           U_state (i,j,bi,bj) =
                0315      &     UtempSt (i,j,bi,bj)
                0316           V_state (i,j,bi,bj) =
                0317      &     VtempSt (i,j,bi,bj)
                0318          ENDDO
                0319         ENDDO
                0320        ENDDO
                0321       ENDDO
                0322
                0323       WRITE(msgBuf,'(A,I5,A)') 'DONE WITH MANUAL CG ADJOINT:',tmpiter,
                0324      & 'ITERS'
                0325        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0326      &                     SQUEEZE_RIGHT , 1)
                0327
                0328 #endif
                0329       RETURN
                0330       END