/[MITgcm]/MITgcm/pkg/seaice/seaice_calc_strainrates.F
ViewVC logotype

Annotation of /MITgcm/pkg/seaice/seaice_calc_strainrates.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.17 - (hide annotations) (download)
Fri Nov 5 08:13:03 2010 UTC (14 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63a, checkpoint63b, checkpoint63c, checkpoint63, checkpoint62o, checkpoint62n, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x
Changes since 1.16: +1 -63 lines
remove SEAICE_OLD_AND_BAD_DISCRETIZATION-code, add a test to
seaice_check that stops the model, when the CPP flag is used.

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

  ViewVC Help
Powered by ViewVC 1.1.22