/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_model.F
ViewVC logotype

Annotation of /MITgcm_contrib/torge/itd/code/seaice_model.F

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


Revision 1.13 - (hide annotations) (download)
Fri May 3 19:21:20 2013 UTC (12 years, 3 months ago) by torge
Branch: MAIN
CVS Tags: HEAD
Changes since 1.12: +1 -143 lines
removing all SEAICE_DEBUG related lines,
which were introduced with SEAICE_ITD development

1 torge 1.13 C $Header: /u/gcmpack/MITgcm_contrib/torge/itd/code/seaice_model.F,v 1.12 2013/05/03 18:59:40 torge Exp $
2 dimitri 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: SEAICE_MODEL
8    
9     C !INTERFACE: ==========================================================
10     SUBROUTINE SEAICE_MODEL( myTime, myIter, myThid )
11    
12     C !DESCRIPTION: \bv
13     C *===========================================================*
14     C | SUBROUTINE SEAICE_MODEL |
15     C | o Time stepping of a dynamic/thermodynamic sea ice model. |
16     C | Dynamics solver: Zhang/Hibler, JGR, 102, 8691-8702, 1997 |
17     C | Thermodynamics: Hibler, MWR, 108, 1943-1973, 1980 |
18     C | Rheology: Hibler, JPO, 9, 815- 846, 1979 |
19     C | Snow: Zhang et al. , JPO, 28, 191- 217, 1998 |
20     C | Parallel forward ice model written by Jinlun Zhang PSC/UW|
21     C | & coupled into MITgcm by Dimitris Menemenlis (JPL) 2/2001|
22     C | zhang@apl.washington.edu / menemenlis@jpl.nasa.gov |
23     C *===========================================================*
24     C *===========================================================*
25     IMPLICIT NONE
26     C \ev
27    
28     C !USES: ===============================================================
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "DYNVARS.h"
32     #include "PARAMS.h"
33     #include "GRID.h"
34     #include "FFIELDS.h"
35     #include "SEAICE_SIZE.h"
36     #include "SEAICE_PARAMS.h"
37     #include "SEAICE.h"
38     #include "SEAICE_TRACER.h"
39     #ifdef ALLOW_EXF
40     # include "EXF_OPTIONS.h"
41     # include "EXF_FIELDS.h"
42     #endif
43     #ifdef ALLOW_AUTODIFF_TAMC
44     # include "tamc.h"
45     #endif
46    
47     C !INPUT PARAMETERS: ===================================================
48     C myTime - Simulation time
49     C myIter - Simulation timestep number
50     C myThid - Thread no. that called this routine.
51     _RL myTime
52     INTEGER myIter
53     INTEGER myThid
54     CEndOfInterface
55    
56     C !LOCAL VARIABLES: ====================================================
57     C i,j,bi,bj :: Loop counters
58 torge 1.9 #ifdef SEAICE_ITD
59     INTEGER IT
60     #endif
61 torge 1.11 #if defined(ALLOW_AUTODIFF_TAMC) || defined(SEAICE_ITD)
62 dimitri 1.1 INTEGER i, j, bi, bj
63     #endif
64     #ifdef ALLOW_SITRACER
65     INTEGER iTr
66     #endif
67     CEOP
68    
69     #ifdef ALLOW_DEBUG
70     IF (debugMode) CALL DEBUG_ENTER( 'SEAICE_MODEL', myThid )
71     #endif
72    
73     C-- Winds are from pkg/exf, which does not update edges.
74     CALL EXCH_UV_AGRID_3D_RL( uwind, vwind, .TRUE., 1, myThid )
75    
76     #ifdef ALLOW_THSICE
77     IF ( useThSice ) THEN
78     C-- Map thSice-variables to HEFF and AREA
79     CALL SEAICE_MAP_THSICE( myTime, myIter, myThid )
80     ENDIF
81     #endif /* ALLOW_THSICE */
82    
83     #ifdef ALLOW_AUTODIFF_TAMC
84     DO bj=myByLo(myThid),myByHi(myThid)
85     DO bi=myBxLo(myThid),myBxHi(myThid)
86     DO j=1-OLy,sNy+OLy
87     DO i=1-OLx,sNx+OLx
88     uIceNm1(i,j,bi,bj) = 0. _d 0
89     vIceNm1(i,j,bi,bj) = 0. _d 0
90     # ifdef ALLOW_SITRACER
91     DO iTr = 1, SItrMaxNum
92     SItrBucket(i,j,bi,bj,iTr) = 0. _d 0
93     ENDDO
94     # endif
95     ENDDO
96     ENDDO
97     ENDDO
98     ENDDO
99     CADJ STORE uwind = comlev1, key=ikey_dynamics, kind=isbyte
100     CADJ STORE vwind = comlev1, key=ikey_dynamics, kind=isbyte
101     CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
102     CADJ STORE heffm = comlev1, key=ikey_dynamics, kind=isbyte
103     CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
104     # ifdef SEAICE_ALLOW_DYNAMICS
105     # ifdef SEAICE_CGRID
106     CADJ STORE seaicemasku = comlev1, key=ikey_dynamics, kind=isbyte
107     CADJ STORE seaicemaskv = comlev1, key=ikey_dynamics, kind=isbyte
108     CADJ STORE fu = comlev1, key=ikey_dynamics, kind=isbyte
109     CADJ STORE fv = comlev1, key=ikey_dynamics, kind=isbyte
110     CADJ STORE uice = comlev1, key=ikey_dynamics, kind=isbyte
111     CADJ STORE vice = comlev1, key=ikey_dynamics, kind=isbyte
112     cphCADJ STORE eta = comlev1, key=ikey_dynamics, kind=isbyte
113     cphCADJ STORE zeta = comlev1, key=ikey_dynamics, kind=isbyte
114     cph(
115     CADJ STORE dwatn = comlev1, key=ikey_dynamics, kind=isbyte
116     cccCADJ STORE press0 = comlev1, key=ikey_dynamics, kind=isbyte
117     cccCADJ STORE taux = comlev1, key=ikey_dynamics, kind=isbyte
118     cccCADJ STORE tauy = comlev1, key=ikey_dynamics, kind=isbyte
119     cccCADJ STORE zmax = comlev1, key=ikey_dynamics, kind=isbyte
120     cccCADJ STORE zmin = comlev1, key=ikey_dynamics, kind=isbyte
121     cph)
122     # ifdef SEAICE_ALLOW_EVP
123     CADJ STORE seaice_sigma1 = comlev1, key=ikey_dynamics, kind=isbyte
124     CADJ STORE seaice_sigma2 = comlev1, key=ikey_dynamics, kind=isbyte
125     CADJ STORE seaice_sigma12 = comlev1, key=ikey_dynamics, kind=isbyte
126     # endif
127     # endif
128     # endif
129     # ifdef ALLOW_SITRACER
130     CADJ STORE siceload = comlev1, key=ikey_dynamics, kind=isbyte
131     CADJ STORE sitracer = comlev1, key=ikey_dynamics, kind=isbyte
132     # endif
133     #endif /* ALLOW_AUTODIFF_TAMC */
134    
135     C solve ice momentum equations and calculate ocean surface stress
136     #ifdef ALLOW_DEBUG
137     IF (debugMode) CALL DEBUG_CALL( 'SEAICE_DYNSOLVER', myThid )
138     #endif
139     #ifdef SEAICE_CGRID
140     CALL TIMER_START('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid)
141     CALL SEAICE_DYNSOLVER ( myTime, myIter, myThid )
142     CALL TIMER_STOP ('SEAICE_DYNSOLVER [SEAICE_MODEL]',myThid)
143     #else
144     CALL TIMER_START('DYNSOLVER [SEAICE_MODEL]',myThid)
145     CALL DYNSOLVER ( myTime, myIter, myThid )
146     CALL TIMER_STOP ('DYNSOLVER [SEAICE_MODEL]',myThid)
147     #endif /* SEAICE_CGRID */
148    
149     C-- Apply ice velocity open boundary conditions
150     #ifdef ALLOW_OBCS
151     # ifndef DISABLE_SEAICE_OBCS
152     IF ( useOBCS ) CALL OBCS_ADJUST_UVICE( uice, vice, myThid )
153     # endif /* DISABLE_SEAICE_OBCS */
154     #endif /* ALLOW_OBCS */
155    
156     #ifdef ALLOW_THSICE
157 torge 1.11 IF ( useThSice ) THEN
158     #ifndef OLD_THSICE_CALL_SEQUENCE
159     #ifdef ALLOW_DEBUG
160     IF (debugMode) CALL DEBUG_CALL( 'THSICE_DO_ADVECT', myThid )
161     #endif
162     CALL THSICE_DO_ADVECT( 0, 0, myTime, myIter, myThid )
163     #endif /* OLD_THSICE_CALL_SEQUENCE */
164     ELSE
165 dimitri 1.1 #endif
166     C-- Only call advection of heff, area, snow, and salt and
167     C-- growth for the generic 0-layer thermodynamics of seaice
168     C-- if useThSice=.false., otherwise the 3-layer Winton thermodynamics
169     C-- (called from DO_OCEANIC_PHYSICS) take care of this
170    
171     C NOW DO ADVECTION and DIFFUSION
172     IF ( SEAICEadvHeff .OR. SEAICEadvArea .OR. SEAICEadvSnow
173     & .OR. SEAICEadvSalt ) THEN
174     #ifdef ALLOW_DEBUG
175     IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ADVDIFF', myThid )
176     #endif
177     CALL SEAICE_ADVDIFF( myTime, myIter, myThid )
178 dimitri 1.2 #ifdef SEAICE_ITD
179     C check that all ice thickness categories meet their limits
180     C (includes Hibler-type ridging)
181     #ifdef ALLOW_DEBUG
182     IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ITD_REDIST', myThid )
183     #endif
184 torge 1.5 DO bj=myByLo(myThid),myByHi(myThid)
185     DO bi=myBxLo(myThid),myBxHi(myThid)
186 torge 1.6 CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
187 torge 1.5 ENDDO
188     ENDDO
189     C update mean ice thickness HEFF and total ice concentration AREA
190     C to match single category values
191     #ifdef ALLOW_DEBUG
192     IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ITD_SUM', myThid )
193     #endif
194     DO bj=myByLo(myThid),myByHi(myThid)
195     DO bi=myBxLo(myThid),myBxHi(myThid)
196 torge 1.6 CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
197 torge 1.5 ENDDO
198     ENDDO
199 dimitri 1.2 #endif
200 dimitri 1.1 ENDIF
201     #ifdef ALLOW_AUTODIFF_TAMC
202     CADJ STORE heffm = comlev1, key=ikey_dynamics, kind=isbyte
203     #endif /* ALLOW_AUTODIFF_TAMC */
204    
205     #ifndef DISABLE_SEAICE_GROWTH
206     C thermodynamics growth
207     C must call growth after calling advection
208     C because of ugly time level business
209     IF ( usePW79thermodynamics ) THEN
210     #ifdef ALLOW_DEBUG
211     IF (debugMode) CALL DEBUG_CALL( 'SEAICE_GROWTH', myThid )
212     #endif
213     CALL SEAICE_GROWTH( myTime, myIter, myThid )
214 dimitri 1.2 #ifdef SEAICE_ITD
215     C redistribute sea ice into proper sea ice category after growth/melt
216     C in case model runs with ice thickness distribution
217     C---+-|--1----+----2----+----3----+----4----+----5----+----6----+----7-|
218     #ifdef ALLOW_DEBUG
219     IF (debugMode) CALL DEBUG_CALL( 'SEAICE_ITD_REDIST', myThid )
220     #endif
221 torge 1.5 DO bj=myByLo(myThid),myByHi(myThid)
222     DO bi=myBxLo(myThid),myBxHi(myThid)
223 torge 1.6 CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
224 torge 1.5 ENDDO
225     ENDDO
226 dimitri 1.2 C store the mean ice thickness in HEFF (for dynamic solver and diagnostics)
227 torge 1.5 DO bj=myByLo(myThid),myByHi(myThid)
228     DO bi=myBxLo(myThid),myBxHi(myThid)
229 torge 1.6 CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
230 torge 1.5 ENDDO
231     ENDDO
232 torge 1.8 #endif
233 dimitri 1.2 C
234 dimitri 1.1 ENDIF
235     #endif /* DISABLE_SEAICE_GROWTH */
236    
237     #ifdef ALLOW_SITRACER
238     # ifdef ALLOW_AUTODIFF_TAMC
239     CADJ STORE sitracer = comlev1, key=ikey_dynamics, kind=isbyte
240     # endif
241     CALL SEAICE_TRACER_PHYS ( myTime, myIter, myThid )
242     #endif
243    
244     C-- Apply ice tracer open boundary conditions
245     #ifdef ALLOW_OBCS
246     # ifndef DISABLE_SEAICE_OBCS
247     IF ( useOBCS ) CALL OBCS_APPLY_SEAICE( myThid )
248     # endif /* DISABLE_SEAICE_OBCS */
249     #endif /* ALLOW_OBCS */
250    
251     C-- Update overlap regions for a bunch of stuff
252     _EXCH_XY_RL( HEFF, myThid )
253     _EXCH_XY_RL( AREA, myThid )
254     _EXCH_XY_RL( HSNOW, myThid )
255     #ifdef SEAICE_VARIABLE_SALINITY
256     _EXCH_XY_RL( HSALT, myThid )
257     #endif
258     #ifdef ALLOW_SITRACER
259     DO iTr = 1, SItrNumInUse
260     _EXCH_XY_RL( SItracer(1-OLx,1-OLy,1,1,iTr),myThid )
261     ENDDO
262     #endif
263     _EXCH_XY_RS(EmPmR, myThid )
264     _EXCH_XY_RS(saltFlux, myThid )
265     _EXCH_XY_RS(Qnet , myThid )
266     #ifdef SHORTWAVE_HEATING
267     _EXCH_XY_RS(Qsw , myThid )
268     #endif /* SHORTWAVE_HEATING */
269     #ifdef ATMOSPHERIC_LOADING
270     IF ( useRealFreshWaterFlux )
271     & _EXCH_XY_RS( sIceLoad, myThid )
272     #endif
273    
274     #ifdef ALLOW_OBCS
275     C-- In case we use scheme with a large stencil that extends into overlap:
276     C no longer needed with the right masking in advection & diffusion S/R.
277     c IF ( useOBCS ) THEN
278     c DO bj=myByLo(myThid),myByHi(myThid)
279     c DO bi=myBxLo(myThid),myBxHi(myThid)
280     c CALL OBCS_COPY_TRACER( HEFF(1-OLx,1-OLy,bi,bj),
281     c I 1, bi, bj, myThid )
282     c CALL OBCS_COPY_TRACER( AREA(1-OLx,1-OLy,bi,bj),
283     c I 1, bi, bj, myThid )
284     c CALL OBCS_COPY_TRACER( HSNOW(1-OLx,1-OLy,bi,bj),
285     c I 1, bi, bj, myThid )
286     #ifdef SEAICE_VARIABLE_SALINITY
287     c CALL OBCS_COPY_TRACER( HSALT(1-OLx,1-OLy,bi,bj),
288     c I 1, bi, bj, myThid )
289     #endif
290     c ENDDO
291     c ENDDO
292     c ENDIF
293     #endif /* ALLOW_OBCS */
294    
295     #ifdef ALLOW_DIAGNOSTICS
296     IF ( useDiagnostics ) THEN
297     C diagnostics for "non-state variables" that are modified by
298     C the seaice model
299     # ifdef ALLOW_EXF
300     CALL DIAGNOSTICS_FILL(UWIND ,'SIuwind ',0,1 ,0,1,1,myThid)
301     CALL DIAGNOSTICS_FILL(VWIND ,'SIvwind ',0,1 ,0,1,1,myThid)
302     # endif
303     CALL DIAGNOSTICS_FILL_RS(FU ,'SIfu ',0,1 ,0,1,1,myThid)
304     CALL DIAGNOSTICS_FILL_RS(FV ,'SIfv ',0,1 ,0,1,1,myThid)
305     CALL DIAGNOSTICS_FILL_RS(EmPmR,'SIempmr ',0,1 ,0,1,1,myThid)
306     CALL DIAGNOSTICS_FILL_RS(Qnet ,'SIqnet ',0,1 ,0,1,1,myThid)
307     CALL DIAGNOSTICS_FILL_RS(Qsw ,'SIqsw ',0,1 ,0,1,1,myThid)
308 torge 1.4 #ifdef SEAICE_ITD
309     CALL DIAGNOSTICS_FILL(HEFFITD ,'SIheffN ',0,nITD,0,1,1,myThid)
310     CALL DIAGNOSTICS_FILL(AREAITD ,'SIareaN ',0,nITD,0,1,1,myThid)
311     #endif
312 dimitri 1.1 ENDIF
313     #endif /* ALLOW_DIAGNOSTICS */
314    
315     #ifdef ALLOW_THSICE
316     C endif .not.useThSice
317     ENDIF
318     #endif /* ALLOW_THSICE */
319     CML This has already been done in seaice_ocean_stress/ostres, so why repeat?
320     CML CALL EXCH_UV_XY_RS(fu,fv,.TRUE.,myThid)
321    
322     #ifdef ALLOW_EXF
323     # ifdef ALLOW_AUTODIFF_TAMC
324     # if (defined (ALLOW_AUTODIFF_MONITOR))
325     CALL EXF_ADJOINT_SNAPSHOTS( 3, myTime, myIter, myThid )
326     # endif
327     # endif
328     #endif
329    
330     #ifdef ALLOW_DEBUG
331     IF (debugMode) CALL DEBUG_LEAVE( 'SEAICE_MODEL', myThid )
332     #endif
333    
334     RETURN
335     END

  ViewVC Help
Powered by ViewVC 1.1.22