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

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

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


Revision 1.1 - (hide annotations) (download)
Fri Apr 28 17:22:44 2017 UTC (8 years, 2 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, HEAD
add sea ice momentum advection code

1 mlosch 1.1 C $Header: $
2     C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5     #ifdef ALLOW_AUTODIFF
6     # include "AUTODIFF_OPTIONS.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: SEAICE_MOM_ADVECTION
11    
12     C !INTERFACE: ==========================================================
13     SUBROUTINE SEAICE_MOM_ADVECTION(
14     I bi,bj,iMin,iMax,jMin,jMax,
15     I uIceLoc, vIceLoc,
16     O gU, gV,
17     I myTime, myIter, myThid )
18     C *==========================================================*
19     C | S/R SEAICE_MOM_ADVECTION |
20     C | o Form the advection of sea ice momentum to be added to |
21     C | the right hand-side of the momentum equation. |
22     C *==========================================================*
23     C | Most of the code is take from S/R MOM_VECINV and reuses |
24     C | code from mom_vecinv and mom_common |
25     C *==========================================================*
26     IMPLICIT NONE
27    
28     C == Global variables ==
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "GRID.h"
33     #include "SEAICE_SIZE.h"
34     #include "SEAICE_PARAMS.h"
35     #ifdef ALLOW_AUTODIFF_TAMC
36     # include "tamc.h"
37     # include "tamc_keys.h"
38     #endif
39    
40     C == Routine arguments ==
41     C bi,bj :: current tile indices
42     C iMin,iMax,jMin,jMax :: loop ranges
43     C uIceLoc ::
44     C vIceLoc ::
45    
46     C gU :: advection tendency (all explicit terms), u component
47     C gV :: advection tendency (all explicit terms), v component
48     C myTime :: current time
49     C myIter :: current time-step number
50     C myThid :: my Thread Id number
51     INTEGER bi,bj
52     INTEGER iMin,iMax,jMin,jMax
53     _RL uIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
54     _RL vIceLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
55     _RL gU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56     _RL gV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57     _RL myTime
58     INTEGER myIter
59     INTEGER myThid
60     CEOP
61    
62     #ifdef SEAICE_ALLOW_MOM_ADVECTION
63    
64     C == Functions ==
65     LOGICAL DIFFERENT_MULTIPLE
66     EXTERNAL DIFFERENT_MULTIPLE
67    
68     C == Local variables ==
69     _RL uCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70     _RL vCf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     _RS hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RS r_hFacZ (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73     _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74     _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75     _RL KE (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76     _RL vort3 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77     C i,j :: Loop counters
78     C k :: surface level index
79     INTEGER i,j,k
80     C later these will be run time parameters
81     CML LOGICAL SEAICEhighOrderVorticity, SEAICEupwindVorticity
82     CML LOGICAL SEAICEuseAbsVorticity,
83     LOGICAL vorticityFlag
84     #ifdef ALLOW_AUTODIFF_TAMC
85     INTEGER imomkey
86     #endif
87    
88     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89    
90     #ifdef ALLOW_AUTODIFF_TAMC
91     act0 = k - 1
92     max0 = Nr
93     act1 = bi - myBxLo(myThid)
94     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
95     act2 = bj - myByLo(myThid)
96     max2 = myByHi(myThid) - myByLo(myThid) + 1
97     act3 = myThid - 1
98     max3 = nTx*nTy
99     act4 = ikey_dynamics - 1
100     imomkey = (act0 + 1)
101     & + act1*max0
102     & + act2*max0*max1
103     & + act3*max0*max1*max2
104     & + act4*max0*max1*max2*max3
105     #endif /* ALLOW_AUTODIFF_TAMC */
106    
107     CML SEAICEselectKEscheme = selectKEscheme
108     CML SEAICEselectVortScheme = selectVortScheme
109     CML SEAICEhighOrderVorticity = highOrderVorticity
110     CML SEAICEupwindVorticity = upwindVorticity
111     CML SEAICEuseAbsVorticity = useAbsVorticity
112     CML SEAICEuseJamartMomAdv = useJamartMomAdv
113    
114     C-- Initialise intermediate terms
115     DO j=1-OLy,sNy+OLy
116     DO i=1-OLx,sNx+OLx
117     uCf(i,j) = 0.
118     vCf(i,j) = 0.
119     gU(i,j) = 0.
120     gV(i,j) = 0.
121     vort3(i,j) = 0.
122     KE(i,j) = 0.
123     #ifdef ALLOW_AUTODIFF
124     hFacZ(i,j) = 0. _d 0
125     #endif
126     ENDDO
127     ENDDO
128    
129     k = 1
130     C-- Calculate open water fraction at vorticity points
131     CALL MOM_CALC_HFACZ(bi,bj,k,hFacZ,r_hFacZ,myThid)
132    
133     C Make local copies of horizontal flow field
134     DO j=1-OLy,sNy+OLy
135     DO i=1-OLx,sNx+OLx
136     uFld(i,j) = uIceLoc(i,j,bi,bj)
137     vFld(i,j) = vIceLoc(i,j,bi,bj)
138     ENDDO
139     ENDDO
140    
141     #ifdef ALLOW_AUTODIFF_TAMC
142     CADJ STORE ufld(:,:) =
143     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
144     CADJ STORE vfld(:,:) =
145     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
146     CADJ STORE hFacZ(:,:) =
147     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
148     CADJ STORE r_hFacZ(:,:) =
149     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
150     #endif
151    
152     CALL MOM_CALC_KE(bi,bj,k,SEAICEselectKEscheme,uFld,vFld,KE,myThid)
153    
154     CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
155    
156     CMLC- calculate absolute vorticity
157     CML IF (useAbsVorticity) THEN
158     CML DO j=1-Oly,sNy+Oly
159     CML DO i=1-Olx,sNx+Olx
160     CML vort3(i,j) = vort3(i,j) + fCoriG(i,j,bi,bj)
161     CML ENDDO
162     CML ENDDO
163     CML ENDIF
164    
165     #ifdef ALLOW_AUTODIFF_TAMC
166     CADJ STORE ke(:,:) =
167     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
168     CADJ STORE vort3(:,:) =
169     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
170     #endif
171    
172     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
173    
174     #ifdef ALLOW_AUTODIFF_TAMC
175     CADJ STORE ucf(:,:) =
176     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
177     CADJ STORE vcf(:,:) =
178     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
179     #endif
180    
181     C-- Horizontal advection of relative (or absolute) vorticity
182     vorticityFlag = SEAICEhighOrderVorticity.OR.SEAICEupwindVorticity
183     IF ( vorticityFlag ) THEN
184     CALL MOM_VI_U_CORIOLIS_C4(bi,bj,k,SEAICEselectVortScheme,
185     & SEAICEhighOrderVorticity,
186     & SEAICEupwindVorticity,
187     & vFld,vort3,r_hFacZ,
188     & uCf,myThid)
189     ELSE
190     CALL MOM_VI_U_CORIOLIS(bi,bj,k,SEAICEselectVortScheme,
191     & SEAICEuseJamartMomAdv,
192     & vFld,vort3,hFacZ,r_hFacZ,
193     & uCf,myThid)
194     ENDIF
195     DO j=jMin,jMax
196     DO i=iMin,iMax
197     gU(i,j) = gU(i,j)+uCf(i,j)
198     ENDDO
199     ENDDO
200     IF ( vorticityFlag ) THEN
201     CALL MOM_VI_V_CORIOLIS_C4(bi,bj,k,SEAICEselectVortScheme,
202     & SEAICEhighOrderVorticity,
203     & SEAICEupwindVorticity,
204     & uFld,vort3,r_hFacZ,
205     & vCf,myThid)
206     ELSE
207     CALL MOM_VI_V_CORIOLIS(bi,bj,k,SEAICEselectVortScheme,
208     & SEAICEuseJamartMomAdv,
209     & uFld,vort3,hFacZ,r_hFacZ,
210     & vCf,myThid)
211     ENDIF
212     DO j=jMin,jMax
213     DO i=iMin,iMax
214     gV(i,j) = gV(i,j)+vCf(i,j)
215     ENDDO
216     ENDDO
217    
218     #ifdef ALLOW_AUTODIFF_TAMC
219     CADJ STORE ucf(:,:) =
220     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
221     CADJ STORE vcf(:,:) =
222     CADJ & comlev1_bibj_lsr, key = imomkey, byte = isbyte
223     #endif
224    
225     #ifdef ALLOW_DIAGNOSTICS
226     IF ( useDiagnostics ) THEN
227     CALL DIAGNOSTICS_FILL(uCf,'SIuAdvZ3',k,1,2,bi,bj,myThid)
228     CALL DIAGNOSTICS_FILL(vCf,'SIvAdvZ3',k,1,2,bi,bj,myThid)
229     ENDIF
230     #endif /* ALLOW_DIAGNOSTICS */
231    
232     C-- Bernoulli term
233     CALL MOM_VI_U_GRAD_KE(bi,bj,k,KE,uCf,myThid)
234     DO j=jMin,jMax
235     DO i=iMin,iMax
236     gU(i,j) = gU(i,j)+uCf(i,j)
237     ENDDO
238     ENDDO
239     CALL MOM_VI_V_GRAD_KE(bi,bj,k,KE,vCf,myThid)
240     DO j=jMin,jMax
241     DO i=iMin,iMax
242     gV(i,j) = gV(i,j)+vCf(i,j)
243     ENDDO
244     ENDDO
245     #ifdef ALLOW_DIAGNOSTICS
246     IF ( useDiagnostics ) THEN
247     CALL DIAGNOSTICS_FILL(uCf,'SIKEx ',k,1,2,bi,bj,myThid)
248     CALL DIAGNOSTICS_FILL(vCf,'SIKEy ',k,1,2,bi,bj,myThid)
249     ENDIF
250     #endif /* ALLOW_DIAGNOSTICS */
251    
252     C-- Set du/dt & dv/dt on boundaries to zero
253     C apply masks for interior (important when we have open boundaries)
254     DO j=jMin,jMax
255     DO i=iMin,iMax
256     gU(i,j) = gU(i,j)*maskInW(i,j,bi,bj)
257     gV(i,j) = gV(i,j)*maskInS(i,j,bi,bj)
258     ENDDO
259     ENDDO
260    
261     #ifdef ALLOW_DIAGNOSTICS
262     IF ( useDiagnostics ) THEN
263     CALL DIAGNOSTICS_FILL(KE, 'SImomKE ',k,1,2,bi,bj,myThid)
264     CALL DIAGNOSTICS_FILL(gU, 'SIuMmAdv',k,1,2,bi,bj,myThid)
265     CALL DIAGNOSTICS_FILL(gV, 'SIvMmAdv',k,1,2,bi,bj,myThid)
266     ENDIF
267     #endif /* ALLOW_DIAGNOSTICS */
268    
269     #endif /* SEAICE_ALLOW_MOM_ADVECTION */
270    
271     RETURN
272     END

  ViewVC Help
Powered by ViewVC 1.1.22