| 1 | mlosch | 1.17 | C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.16 2010/03/16 19:21:31 gforget Exp $ | 
| 2 | mlosch | 1.1 | C $Name:  $ | 
| 3 |  |  |  | 
| 4 |  |  | #include "SEAICE_OPTIONS.h" | 
| 5 |  |  |  | 
| 6 |  |  | CStartOfInterface | 
| 7 | jmc | 1.8 | SUBROUTINE SEAICE_CALC_STRAINRATES( | 
| 8 | mlosch | 1.1 | I     uFld, vFld, | 
| 9 | mlosch | 1.12 | O     e11Loc, e22Loc, e12Loc, | 
| 10 | mlosch | 1.14 | I     iStep, myTime, myIter, myThid ) | 
| 11 | mlosch | 1.1 | C     /==========================================================\ | 
| 12 |  |  | C     | SUBROUTINE  SEAICE_CALC_STRAINRATES                      | | 
| 13 |  |  | C     | o compute strain rates from ice velocities               | | 
| 14 |  |  | C     |==========================================================| | 
| 15 |  |  | C     | written by Martin Losch, Apr 2007                        | | 
| 16 |  |  | C     \==========================================================/ | 
| 17 |  |  | IMPLICIT NONE | 
| 18 |  |  |  | 
| 19 |  |  | C     === Global variables === | 
| 20 |  |  | #include "SIZE.h" | 
| 21 |  |  | #include "EEPARAMS.h" | 
| 22 |  |  | #include "PARAMS.h" | 
| 23 |  |  | #include "GRID.h" | 
| 24 |  |  | #include "SEAICE_PARAMS.h" | 
| 25 | mlosch | 1.11 | #include "SEAICE.h" | 
| 26 | mlosch | 1.1 |  | 
| 27 |  |  | #ifdef ALLOW_AUTODIFF_TAMC | 
| 28 |  |  | # include "tamc.h" | 
| 29 |  |  | #endif | 
| 30 |  |  |  | 
| 31 |  |  | C     === Routine arguments === | 
| 32 | jmc | 1.8 | C     iStep  :: Sub-time-step number | 
| 33 |  |  | C     myTime :: Simulation time | 
| 34 |  |  | C     myIter :: Simulation timestep number | 
| 35 |  |  | C     myThid :: My Thread Id. number | 
| 36 |  |  | INTEGER iStep | 
| 37 |  |  | _RL     myTime | 
| 38 |  |  | INTEGER myIter | 
| 39 | mlosch | 1.1 | INTEGER myThid | 
| 40 |  |  | C     ice velocities | 
| 41 | mlosch | 1.14 | _RL uFld   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) | 
| 42 |  |  | _RL vFld   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy) | 
| 43 | mlosch | 1.1 | C     strain rate tensor | 
| 44 | mlosch | 1.12 | _RL e11Loc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) | 
| 45 |  |  | _RL e22Loc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) | 
| 46 |  |  | _RL e12Loc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) | 
| 47 | mlosch | 1.1 | CEndOfInterface | 
| 48 |  |  |  | 
| 49 |  |  | #ifdef SEAICE_CGRID | 
| 50 |  |  | #ifdef SEAICE_ALLOW_DYNAMICS | 
| 51 |  |  | C     === Local variables === | 
| 52 |  |  | C     i,j,bi,bj - Loop counters | 
| 53 |  |  | INTEGER i, j, bi, bj | 
| 54 | mlosch | 1.15 | C     hFacU, hFacV - determine the no-slip boundary condition | 
| 55 | mlosch | 1.2 | INTEGER k | 
| 56 | mlosch | 1.11 | _RS hFacU, hFacV, noSlipFac | 
| 57 | mlosch | 1.15 | C     auxillary variables that help writing code that | 
| 58 |  |  | C     vectorizes even after TAFization | 
| 59 |  |  | _RL dudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 60 |  |  | _RL dvdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 61 |  |  | _RL dudy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 62 |  |  | _RL dvdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 63 |  |  | _RL uave (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 64 |  |  | _RL vave (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 65 | mlosch | 1.2 |  | 
| 66 | mlosch | 1.4 | k = 1 | 
| 67 | mlosch | 1.11 | noSlipFac = 0. _d 0 | 
| 68 |  |  | IF ( SEAICE_no_slip ) noSlipFac = 1. _d 0 | 
| 69 | mlosch | 1.1 | C | 
| 70 | mlosch | 1.11 | DO bj=myByLo(myThid),myByHi(myThid) | 
| 71 |  |  | DO bi=myBxLo(myThid),myBxHi(myThid) | 
| 72 | mlosch | 1.15 | C     abbreviations on C-points, need to do them in separate loops | 
| 73 |  |  | C     for vectorization | 
| 74 | mlosch | 1.11 | DO j=1-Oly,sNy+Oly-1 | 
| 75 |  |  | DO i=1-Olx,sNx+Olx-1 | 
| 76 | mlosch | 1.15 | dudx(I,J) = _recip_dxF(I,J,bi,bj) * | 
| 77 | mlosch | 1.14 | &         (uFld(I+1,J,bi,bj)-uFld(I,J,bi,bj)) | 
| 78 | mlosch | 1.15 | uave(I,J) = 0.5 _d 0 * (uFld(I,J,bi,bj)+uFld(I+1,J,bi,bj)) | 
| 79 |  |  | ENDDO | 
| 80 |  |  | ENDDO | 
| 81 |  |  | DO j=1-Oly,sNy+Oly-1 | 
| 82 |  |  | DO i=1-Olx,sNx+Olx-1 | 
| 83 |  |  | dvdy(I,J) = _recip_dyF(I,J,bi,bj) * | 
| 84 | mlosch | 1.14 | &         (vFld(I,J+1,bi,bj)-vFld(I,J,bi,bj)) | 
| 85 | mlosch | 1.15 | vave(I,J) = 0.5 _d 0 * (vFld(I,J,bi,bj)+vFld(I,J+1,bi,bj)) | 
| 86 |  |  | ENDDO | 
| 87 |  |  | ENDDO | 
| 88 |  |  | C     evaluate strain rates at C-points | 
| 89 |  |  | DO j=1-Oly,sNy+Oly-1 | 
| 90 |  |  | DO i=1-Olx,sNx+Olx-1 | 
| 91 |  |  | e11Loc(I,J,bi,bj) = dudx(I,J) + vave(I,J) * k2AtC(I,J,bi,bj) | 
| 92 |  |  | e22Loc(I,J,bi,bj) = dvdy(I,J) + uave(I,J) * k1AtC(I,J,bi,bj) | 
| 93 | mlosch | 1.11 | ENDDO | 
| 94 |  |  | ENDDO | 
| 95 | mlosch | 1.15 | C     abbreviations at Z-points, need to do them in separate loops | 
| 96 |  |  | C     for vectorization | 
| 97 |  |  | DO j=1-Oly+1,sNy+Oly | 
| 98 |  |  | DO i=1-Olx+1,sNx+Olx | 
| 99 |  |  | dudy(I,J) = ( uFld(I,J,bi,bj) - uFld(I  ,J-1,bi,bj) ) | 
| 100 |  |  | &         * _recip_dyU(I,J,bi,bj) | 
| 101 |  |  | uave(I,J) = 0.5 _d 0 * (uFld(I,J,bi,bj)+uFld(I  ,J-1,bi,bj)) | 
| 102 |  |  | ENDDO | 
| 103 |  |  | ENDDO | 
| 104 |  |  | DO j=1-Oly+1,sNy+Oly | 
| 105 |  |  | DO i=1-Olx+1,sNx+Olx | 
| 106 |  |  | dvdx(I,J) = ( vFld(I,J,bi,bj) - vFld(I-1,J  ,bi,bj) ) | 
| 107 |  |  | &         * _recip_dxV(I,J,bi,bj) | 
| 108 |  |  | vave(I,J) = 0.5 _d 0 * (vFld(I,J,bi,bj)+vFld(I-1,J  ,bi,bj)) | 
| 109 |  |  | ENDDO | 
| 110 |  |  | ENDDO | 
| 111 |  |  | C     evaluate strain rates at Z-points | 
| 112 | mlosch | 1.11 | DO j=1-Oly+1,sNy+Oly | 
| 113 |  |  | DO i=1-Olx+1,sNx+Olx | 
| 114 |  |  | hFacU = _maskW(i,j,k,bi,bj) - _maskW(i,j-1,k,bi,bj) | 
| 115 |  |  | hFacV = _maskS(i,j,k,bi,bj) - _maskS(i-1,j,k,bi,bj) | 
| 116 | mlosch | 1.15 | e12Loc(I,J,bi,bj) = 0.5 _d 0 * ( | 
| 117 |  |  | &         dudy(I,J) + dvdx(I,J) | 
| 118 |  |  | &         - k1AtZ(I,J,bi,bj) * vave(I,J) | 
| 119 |  |  | &         - k2AtZ(I,J,bi,bj) * uave(I,J) | 
| 120 | mlosch | 1.11 | &         ) | 
| 121 |  |  | &         *maskC(I  ,J  ,k,bi,bj)*maskC(I-1,J  ,k,bi,bj) | 
| 122 |  |  | &         *maskC(I  ,J-1,k,bi,bj)*maskC(I-1,J-1,k,bi,bj) | 
| 123 |  |  | &         + 2.0 _d 0 * noSlipFac * ( | 
| 124 | mlosch | 1.15 | &           2.0 _d 0 * uave(I,J) * _recip_dyU(I,J,bi,bj) * hFacU | 
| 125 |  |  | &         + 2.0 _d 0 * vave(I,J) * _recip_dxV(I,J,bi,bj) * hFacV | 
| 126 | mlosch | 1.11 | &         ) | 
| 127 |  |  | C     no slip at the boundary implies u(j)+u(j-1)=0 and v(i)+v(i-1)=0 | 
| 128 |  |  | C     accross the boundary; this is already accomplished by masking so | 
| 129 |  |  | C     that the following lines are not necessary | 
| 130 | mlosch | 1.15 | c$$$     &         - hFacV * k1AtZ(I,J,bi,bj) * vave(I,J) | 
| 131 |  |  | c$$$     &         - hFacU * k2AtZ(I,J,bi,bj) * uave(I,J) | 
| 132 | mlosch | 1.11 | ENDDO | 
| 133 |  |  | ENDDO | 
| 134 |  |  |  | 
| 135 |  |  | ENDDO | 
| 136 |  |  | ENDDO | 
| 137 | gforget | 1.16 |  | 
| 138 |  |  | #ifdef ALLOW_AUTODIFF_TAMC | 
| 139 |  |  | #ifdef SEAICE_DYN_STABLE_ADJOINT | 
| 140 |  |  | cgf zero out adjoint fields to stabilize pkg/seaice dyna. adjoint | 
| 141 |  |  | CALL ZERO_ADJ( 1, e11Loc, myThid) | 
| 142 |  |  | CALL ZERO_ADJ( 1, e12Loc, myThid) | 
| 143 |  |  | CALL ZERO_ADJ( 1, e22Loc, myThid) | 
| 144 |  |  | #endif | 
| 145 |  |  | #endif /* ALLOW_AUTODIFF_TAMC */ | 
| 146 |  |  |  | 
| 147 | mlosch | 1.1 | #endif /* SEAICE_ALLOW_DYNAMICS */ | 
| 148 |  |  | #endif /* SEAICE_CGRID */ | 
| 149 |  |  | RETURN | 
| 150 |  |  | END |