/[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.21 - (hide annotations) (download)
Mon Oct 20 03:20:57 2014 UTC (10 years, 9 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65g
Changes since 1.20: +4 -1 lines 
- ECCO_OPTIONS.h is needed when including ecco_cost.h, ecco.h
- AUTODIFF_OPTIONS.h is needed when including tamc.h, tamc_keys.h
- CTRL_OPTIONS.h is needed when including ctrl.h, etc

- pkg/seaice/seaice_cost*.F : clean up CPP brackets
- SEAICE_SIZE.h : replace ALLOW_AUTODIFF_TAMC with ALLOW_AUTODIFF to
  avoid needing AUTODIFF_OPTIONS.h anytime SEAICE_SIZE.h is included
  (it seems that THSICE_SIZE.h, PTRACERS_SIZE.h have the same issue...)
1 gforget 1.21 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_calc_strainrates.F,v 1.20 2013/02/28 17:25:41 mlosch Exp $
2 mlosch 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5 jmc 1.18 #ifdef ALLOW_OBCS
6     # include "OBCS_OPTIONS.h"
7     #else
8     # define OBCS_UVICE_OLD
9     #endif
10 gforget 1.21 #ifdef ALLOW_AUTODIFF
11     # include "AUTODIFF_OPTIONS.h"
12     #endif
13 mlosch 1.1
14 jmc 1.18 CBOP
15     C !ROUTINE: SEAICE_CALC_STRAINRATES
16     C !INTERFACE:
17 jmc 1.8 SUBROUTINE SEAICE_CALC_STRAINRATES(
18 mlosch 1.1 I uFld, vFld,
19 mlosch 1.12 O e11Loc, e22Loc, e12Loc,
20 mlosch 1.14 I iStep, myTime, myIter, myThid )
21 jmc 1.18
22     C !DESCRIPTION: \bv
23     C *==========================================================*
24     C | SUBROUTINE SEAICE_CALC_STRAINRATES
25     C | o compute strain rates from ice velocities
26     C *==========================================================*
27     C | written by Martin Losch, Apr 2007
28     C *==========================================================*
29     C \ev
30    
31     C !USES:
32 mlosch 1.1 IMPLICIT NONE
33    
34     C === Global variables ===
35     #include "SIZE.h"
36     #include "EEPARAMS.h"
37     #include "PARAMS.h"
38     #include "GRID.h"
39 jmc 1.19 #include "SEAICE_SIZE.h"
40 mlosch 1.1 #include "SEAICE_PARAMS.h"
41 mlosch 1.11 #include "SEAICE.h"
42 mlosch 1.1
43     #ifdef ALLOW_AUTODIFF_TAMC
44     # include "tamc.h"
45     #endif
46    
47 jmc 1.18 C !INPUT/OUTPUT PARAMETERS:
48 mlosch 1.1 C === Routine arguments ===
49 jmc 1.18 C uFld :: ice velocity, u-component
50     C vFld :: ice velocity, v-component
51     C e11Loc :: strain rate tensor, component 1,1
52     C e22Loc :: strain rate tensor, component 2,2
53     C e12Loc :: strain rate tensor, component 1,2
54 jmc 1.8 C iStep :: Sub-time-step number
55     C myTime :: Simulation time
56     C myIter :: Simulation timestep number
57     C myThid :: My Thread Id. number
58 jmc 1.19 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
59     _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
60 mlosch 1.12 _RL e11Loc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
61     _RL e22Loc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
62     _RL e12Loc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
63 jmc 1.18 INTEGER iStep
64     _RL myTime
65     INTEGER myIter
66     INTEGER myThid
67     CEOP
68 mlosch 1.1
69     #ifdef SEAICE_CGRID
70     #ifdef SEAICE_ALLOW_DYNAMICS
71 jmc 1.18 C !LOCAL VARIABLES:
72 mlosch 1.1 C === Local variables ===
73 jmc 1.18 C i,j,bi,bj :: Loop counters
74 mlosch 1.1 INTEGER i, j, bi, bj
75 jmc 1.18 C hFacU, hFacV :: determine the no-slip boundary condition
76 mlosch 1.2 INTEGER k
77 mlosch 1.11 _RS hFacU, hFacV, noSlipFac
78 mlosch 1.15 C auxillary variables that help writing code that
79     C vectorizes even after TAFization
80     _RL dudx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81     _RL dvdy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82     _RL dudy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83     _RL dvdx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84     _RL uave (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85     _RL vave (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86 mlosch 1.2
87 mlosch 1.4 k = 1
88 mlosch 1.11 noSlipFac = 0. _d 0
89     IF ( SEAICE_no_slip ) noSlipFac = 1. _d 0
90 mlosch 1.20 C in order repoduce results before fixing a bug in r1.20 comment out
91     C the following line
92     CML IF ( SEAICE_no_slip ) noSlipFac = 2. _d 0
93 mlosch 1.1 C
94 mlosch 1.11 DO bj=myByLo(myThid),myByHi(myThid)
95     DO bi=myBxLo(myThid),myBxHi(myThid)
96 mlosch 1.15 C abbreviations on C-points, need to do them in separate loops
97     C for vectorization
98 jmc 1.19 DO j=1-OLy,sNy+OLy-1
99     DO i=1-OLx,sNx+OLx-1
100 jmc 1.18 dudx(i,j) = _recip_dxF(i,j,bi,bj) *
101     & (uFld(i+1,j,bi,bj)-uFld(i,j,bi,bj))
102     uave(i,j) = 0.5 _d 0 * (uFld(i,j,bi,bj)+uFld(i+1,j,bi,bj))
103 mlosch 1.15 ENDDO
104     ENDDO
105 jmc 1.19 DO j=1-OLy,sNy+OLy-1
106     DO i=1-OLx,sNx+OLx-1
107 jmc 1.18 dvdy(i,j) = _recip_dyF(i,j,bi,bj) *
108     & (vFld(i,j+1,bi,bj)-vFld(i,j,bi,bj))
109     vave(i,j) = 0.5 _d 0 * (vFld(i,j,bi,bj)+vFld(i,j+1,bi,bj))
110 mlosch 1.15 ENDDO
111     ENDDO
112     C evaluate strain rates at C-points
113 jmc 1.19 DO j=1-OLy,sNy+OLy-1
114     DO i=1-OLx,sNx+OLx-1
115 jmc 1.18 e11Loc(i,j,bi,bj) = dudx(i,j) + vave(i,j) * k2AtC(i,j,bi,bj)
116     e22Loc(i,j,bi,bj) = dvdy(i,j) + uave(i,j) * k1AtC(i,j,bi,bj)
117     ENDDO
118     ENDDO
119     #ifndef OBCS_UVICE_OLD
120     C-- for OBCS: assume no gradient beyong OB
121 jmc 1.19 DO j=1-OLy,sNy+OLy-1
122     DO i=1-OLx,sNx+OLx-1
123 jmc 1.18 e11Loc(i,j,bi,bj) = e11Loc(i,j,bi,bj)*maskInC(i,j,bi,bj)
124     e22Loc(i,j,bi,bj) = e22Loc(i,j,bi,bj)*maskInC(i,j,bi,bj)
125 mlosch 1.11 ENDDO
126     ENDDO
127 jmc 1.18 #endif /* OBCS_UVICE_OLD */
128    
129 mlosch 1.15 C abbreviations at Z-points, need to do them in separate loops
130     C for vectorization
131 jmc 1.19 DO j=1-OLy+1,sNy+OLy
132     DO i=1-OLx+1,sNx+OLx
133 jmc 1.18 dudy(i,j) = ( uFld(i,j,bi,bj) - uFld(i ,j-1,bi,bj) )
134     & * _recip_dyU(i,j,bi,bj)
135     uave(i,j) = 0.5 _d 0 * (uFld(i,j,bi,bj)+uFld(i ,j-1,bi,bj))
136 mlosch 1.15 ENDDO
137     ENDDO
138 jmc 1.19 DO j=1-OLy+1,sNy+OLy
139     DO i=1-OLx+1,sNx+OLx
140 jmc 1.18 dvdx(i,j) = ( vFld(i,j,bi,bj) - vFld(i-1,j ,bi,bj) )
141     & * _recip_dxV(i,j,bi,bj)
142     vave(i,j) = 0.5 _d 0 * (vFld(i,j,bi,bj)+vFld(i-1,j ,bi,bj))
143 mlosch 1.15 ENDDO
144     ENDDO
145     C evaluate strain rates at Z-points
146 jmc 1.19 DO j=1-OLy+1,sNy+OLy
147     DO i=1-OLx+1,sNx+OLx
148 mlosch 1.11 hFacU = _maskW(i,j,k,bi,bj) - _maskW(i,j-1,k,bi,bj)
149     hFacV = _maskS(i,j,k,bi,bj) - _maskS(i-1,j,k,bi,bj)
150 jmc 1.18 e12Loc(i,j,bi,bj) = 0.5 _d 0 * (
151     & dudy(i,j) + dvdx(i,j)
152     & - k1AtZ(i,j,bi,bj) * vave(i,j)
153     & - k2AtZ(i,j,bi,bj) * uave(i,j)
154 mlosch 1.11 & )
155 jmc 1.18 & *maskC(i ,j ,k,bi,bj)*maskC(i-1,j ,k,bi,bj)
156     & *maskC(i ,j-1,k,bi,bj)*maskC(i-1,j-1,k,bi,bj)
157 mlosch 1.20 & + noSlipFac * (
158 jmc 1.18 & 2.0 _d 0 * uave(i,j) * _recip_dyU(i,j,bi,bj) * hFacU
159     & + 2.0 _d 0 * vave(i,j) * _recip_dxV(i,j,bi,bj) * hFacV
160 mlosch 1.11 & )
161     C no slip at the boundary implies u(j)+u(j-1)=0 and v(i)+v(i-1)=0
162     C accross the boundary; this is already accomplished by masking so
163     C that the following lines are not necessary
164 jmc 1.18 c$$$ & - hFacV * k1AtZ(i,j,bi,bj) * vave(i,j)
165     c$$$ & - hFacU * k2AtZ(i,j,bi,bj) * uave(i,j)
166 mlosch 1.11 ENDDO
167     ENDDO
168    
169     ENDDO
170     ENDDO
171 gforget 1.16
172     #ifdef ALLOW_AUTODIFF_TAMC
173     #ifdef SEAICE_DYN_STABLE_ADJOINT
174     cgf zero out adjoint fields to stabilize pkg/seaice dyna. adjoint
175     CALL ZERO_ADJ( 1, e11Loc, myThid)
176     CALL ZERO_ADJ( 1, e12Loc, myThid)
177     CALL ZERO_ADJ( 1, e22Loc, myThid)
178     #endif
179     #endif /* ALLOW_AUTODIFF_TAMC */
180    
181 mlosch 1.1 #endif /* SEAICE_ALLOW_DYNAMICS */
182     #endif /* SEAICE_CGRID */
183     RETURN
184     END
  ViewVC Help
Powered by ViewVC 1.1.22