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

Contents of /MITgcm/pkg/seaice/seaice_reg_ridge.F

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


Revision 1.4 - (show annotations) (download)
Tue Jul 8 18:56:29 2014 UTC (9 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64z, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65d, checkpoint65e, checkpoint65
Changes since 1.3: +12 -13 lines
remove unused variables

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_reg_ridge.F,v 1.3 2014/06/10 07:17:05 mlosch Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5 #ifdef ALLOW_EXF
6 # include "EXF_OPTIONS.h"
7 #endif
8
9 CBOP
10 C !ROUTINE: SEAICE_REG_RIDGE
11 C !INTERFACE:
12 SUBROUTINE SEAICE_REG_RIDGE( myTime, myIter, myThid )
13 C !DESCRIPTION: \bv
14 C *=================================================================*
15 C | SUBROUTINE seaice_reg_ridge
16 C | o this routine has two purposes:
17 C | (1) clean up after advection (undershoots etc.);
18 C | after advection, the sea ice variables may have unphysical
19 C | values, e.g. < 0 or very thin ice, that are regularized
20 C | here.
21 C | (2) driver for ice ridging;
22 C | concentration as a special case may be > 1 in convergent
23 C | motion and a ridging algorithm redistributes the ice to
24 C | limit the concentration to 1.
25 C | o called after S/R seaice_advdiff
26 C *=================================================================*
27 C \ev
28
29 C !USES:
30 IMPLICIT NONE
31 C === Global variables ===
32 #include "SIZE.h"
33 #include "EEPARAMS.h"
34 #include "PARAMS.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_FIELDS.h"
41 #endif
42 #ifdef ALLOW_AUTODIFF_TAMC
43 # include "tamc.h"
44 #endif /* ALLOW_AUTODIFF_TAMC */
45
46 C !INPUT/OUTPUT PARAMETERS:
47 C === Routine arguments ===
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, myThid
53
54 #ifdef ALLOW_SEAICE
55 C !LOCAL VARIABLES:
56 C === Local variables ===
57 C i,j,bi,bj :: Loop counters
58 INTEGER i, j, bi, bj
59 C number of surface interface layer
60 C IT :: ice thickness category index (ITD and SEAICE_multDim code)
61 INTEGER IT
62 C reciprocal of time step
63 _RL recip_deltaTtherm
64 C temporary variables available for the various computations
65 _RL tmpscal1, tmpscal2
66 #ifdef SEAICE_ITD
67 _RL tmpscal1itd(1:sNx,1:sNy), tmpscal2itd(1:sNx,1:sNy)
68 _RL tmpscal3itd(1:sNx,1:sNy)
69 C reciprocal number of ice classes nITD
70 _RL recip_nitd
71 #endif /* SEAICE_ITD */
72 #ifdef ALLOW_DIAGNOSTICS
73 C Helper variables for diagnostics
74 _RL DIAGarrayA (1:sNx,1:sNy)
75 #endif /* ALLOW_DIAGNOSTICS */
76 CEOP
77
78 C
79 C === Routine body ===
80 C
81 recip_deltaTtherm = ONE / SEAICE_deltaTtherm
82
83 DO bj=myByLo(myThid),myByHi(myThid)
84 DO bi=myBxLo(myThid),myBxHi(myThid)
85
86 #ifdef ALLOW_AUTODIFF_TAMC
87 act1 = bi - myBxLo(myThid)
88 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
89 act2 = bj - myByLo(myThid)
90 max2 = myByHi(myThid) - myByLo(myThid) + 1
91 act3 = myThid - 1
92 max3 = nTx*nTy
93 act4 = ikey_dynamics - 1
94 iicekey = (act1 + 1) + act2*max1
95 & + act3*max1*max2
96 & + act4*max1*max2*max3
97 #endif /* ALLOW_AUTODIFF_TAMC */
98
99 DO J=1-OLy,sNy+OLy
100 DO I=1-OLx,sNx+OLx
101 d_HEFFbyNEG(I,J,bi,bj) = 0.0 _d 0
102 d_HSNWbyNEG(I,J,bi,bj) = 0.0 _d 0
103 #ifdef EXF_SEAICE_FRACTION
104 d_AREAbyRLX(I,J,bi,bj) = 0.0 _d 0
105 d_HEFFbyRLX(I,J,bi,bj) = 0.0 _d 0
106 #endif /* EXF_SEAICE_FRACTION */
107 #ifdef SEAICE_VARIABLE_SALINITY
108 saltFluxAdjust(I,J,bi,bj) = 0.0 _d 0
109 #endif /* SEAICE_VARIABLE_SALINITY */
110 ENDDO
111 ENDDO
112
113 C =====================================================================
114 C ========== PART 1: treat pathological cases (post advdiff) ==========
115 C =====================================================================
116
117 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
118 Cgf no dependency through pathological cases treatment
119 IF ( SEAICEadjMODE.EQ.0 ) THEN
120 #endif
121
122 #ifdef EXF_SEAICE_FRACTION
123 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
124 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
125 C-- (0) relax sea ice concentration towards observation
126 IF ( SEAICE_tauAreaObsRelax .GT. zeroRL ) THEN
127 DO J=1,sNy
128 DO I=1,sNx
129 IF ( exf_iceFraction(I,J,bi,bj).GT.AREA(I,J,bi,bj) ) THEN
130 d_AREAbyRLX(i,j,bi,bj) =
131 & SEAICE_deltaTtherm/SEAICE_tauAreaObsRelax
132 & * (exf_iceFraction(I,J,bi,bj) - AREA(I,J,bi,bj))
133 ENDIF
134 IF ( exf_iceFraction(I,J,bi,bj).GT.zeroRS .AND.
135 & AREA(I,J,bi,bj).EQ.0. _d 0) THEN
136 C d_HEFFbyRLX(i,j,bi,bj) = 1. _d 1 * siEps * d_AREAbyRLX(i,j,bi,bj)
137 d_HEFFbyRLX(i,j,bi,bj) = 1. _d 1 * siEps
138 ENDIF
139 #ifdef SEAICE_ITD
140 AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
141 & + d_AREAbyRLX(i,j,bi,bj)
142 HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
143 & + d_HEFFbyRLX(i,j,bi,bj)
144 #endif /* SEAICE_ITD */
145 AREA(I,J,bi,bj) = AREA(I,J,bi,bj) + d_AREAbyRLX(i,j,bi,bj)
146 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + d_HEFFbyRLX(i,j,bi,bj)
147 ENDDO
148 ENDDO
149 ENDIF
150 #endif /* EXF_SEAICE_FRACTION */
151
152 C-- (1) treat the case of negative values:
153
154 #ifdef SEAICE_ITD
155 DO IT=1,SEAICE_multDim
156 DO J=1,sNy
157 DO I=1,sNx
158 tmpscal1=0. _d 0
159 tmpscal2=0. _d 0
160 tmpscal1=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
161 HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal1
162 d_HEFFbyNEG(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)+tmpscal1
163 tmpscal2=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
164 HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal2
165 d_HSNWbyNEG(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)+tmpscal2
166 AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
167 C AREA, HEFF, and HSNOW will be updated at end of PART 1
168 C by calling SEAICE_ITD_SUM
169 ENDDO
170 ENDDO
171 ENDDO
172 C update mean thicknesses HEFF and HSNOW and total ice
173 C concentration AREA to match single category values
174 CALL SEAICE_ITD_SUM ( bi, bj, myTime, myIter, myThid )
175 #else /* ndef SEAICE_ITD */
176 #ifdef ALLOW_AUTODIFF_TAMC
177 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
178 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
179 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
180 #endif /* ALLOW_AUTODIFF_TAMC */
181 DO J=1,sNy
182 DO I=1,sNx
183 d_HEFFbyNEG(I,J,bi,bj)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
184 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J,bi,bj)
185 d_HSNWbyNEG(I,J,bi,bj)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
186 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J,bi,bj)
187 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
188 ENDDO
189 ENDDO
190 #endif /* SEAICE_ITD */
191
192 C-- (2) treat the case of very thin ice:
193
194 #ifdef SEAICE_ITD
195 C Here we risk that even though HEFF may be larger than siEps (=1e-5)
196 C HEFFITD can have classes with very small (< siEps) non-zero ice volume.
197 C We avoid applying the correction to each class because that leads to
198 C funny structures in the net heat and freshwater flux into the ocean.
199 C Let us keep our fingers crossed, that the model will be benign!
200 DO IT=1,SEAICE_multDim
201 DO J=1,sNy
202 DO I=1,sNx
203 IF (HEFF(I,J,bi,bj).LE.siEps) THEN
204 HEFFITD(I,J,IT,bi,bj) = 0. _d 0
205 HSNOWITD(I,J,IT,bi,bj) = 0. _d 0
206 ENDIF
207 ENDDO
208 ENDDO
209 ENDDO
210 #endif /* SEAICE_ITD */
211 #ifdef ALLOW_AUTODIFF_TAMC
212 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
213 #endif /* ALLOW_AUTODIFF_TAMC */
214 DO J=1,sNy
215 DO I=1,sNx
216 tmpscal1=0. _d 0
217 tmpscal2=0. _d 0
218 IF (HEFF(I,J,bi,bj).LE.siEps) THEN
219 tmpscal1=-HEFF(I,J,bi,bj)
220 tmpscal2=-HSNOW(I,J,bi,bj)
221 DO IT=1,SEAICE_multDim
222 TICES(I,J,IT,bi,bj)=celsius2K
223 ENDDO
224 ENDIF
225 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal1
226 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal2
227 d_HEFFbyNEG(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)+tmpscal1
228 d_HSNWbyNEG(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)+tmpscal2
229 ENDDO
230 ENDDO
231
232 C-- (3) treat the case of area but no ice/snow:
233
234 #ifdef SEAICE_ITD
235 DO IT=1,SEAICE_multDim
236 DO J=1,sNy
237 DO I=1,sNx
238 IF ( (HEFFITD(I,J,IT,bi,bj) .EQ.0. _d 0).AND.
239 & (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
240 & AREAITD(I,J,IT,bi,bj)=0. _d 0
241 ENDDO
242 ENDDO
243 ENDDO
244 #else /* ndef SEAICE_ITD */
245 #ifdef ALLOW_AUTODIFF_TAMC
246 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
247 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
248 #endif /* ALLOW_AUTODIFF_TAMC */
249 DO J=1,sNy
250 DO I=1,sNx
251 IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
252 & (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
253 ENDDO
254 ENDDO
255 #endif /* SEAICE_ITD */
256
257 C-- (4) treat the case of very small area:
258
259 #ifndef DISABLE_AREA_FLOOR
260 #ifdef SEAICE_ITD
261 recip_nitd = 1. _d 0 / float(SEAICE_multDim)
262 DO IT=1,SEAICE_multDim
263 DO J=1,sNy
264 DO I=1,sNx
265 IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
266 & (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
267 C SEAICE_area_floor*SEAICE_multDim cannot be allowed to exceed 1
268 C hence use SEAICE_area_floor devided by SEAICE_multDim
269 C (or install a warning in e.g. seaice_readparms.F)
270 AREAITD(I,J,IT,bi,bj)=
271 & MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor*recip_nitd)
272 ENDIF
273 ENDDO
274 ENDDO
275 ENDDO
276 #else /* ndef SEAICE_ITD */
277 #ifdef ALLOW_AUTODIFF_TAMC
278 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
279 #endif /* ALLOW_AUTODIFF_TAMC */
280 DO J=1,sNy
281 DO I=1,sNx
282 IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
283 AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
284 ENDIF
285 ENDDO
286 ENDDO
287 #endif /* SEAICE_ITD */
288 #endif /* DISABLE_AREA_FLOOR */
289
290 C (5) treat sea ice salinity pathological cases
291 #ifdef SEAICE_VARIABLE_SALINITY
292 #ifdef ALLOW_AUTODIFF_TAMC
293 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
294 CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
295 #endif /* ALLOW_AUTODIFF_TAMC */
296 DO J=1,sNy
297 DO I=1,sNx
298 IF ( (HSALT(I,J,bi,bj) .LT. 0.0).OR.
299 & (HEFF(I,J,bi,bj) .EQ. 0.0) ) THEN
300 saltFluxAdjust(I,J,bi,bj) = - HEFFM(I,J,bi,bj) *
301 & HSALT(I,J,bi,bj) * recip_deltaTtherm
302 HSALT(I,J,bi,bj) = 0.0 _d 0
303 ENDIF
304 ENDDO
305 ENDDO
306 #endif /* SEAICE_VARIABLE_SALINITY */
307
308 C =====================================================================
309 C ========== PART 2: ridging algorithm ===============================
310 C =====================================================================
311
312 C treat case of excessive ice cover, e.g., due to ridging:
313
314 #ifdef SEAICE_ITD
315
316 C catch up with item (2) that involves category sums AREA and HEFF
317 DO J=1,sNy
318 DO I=1,sNx
319 tmpscal1itd(i,j) = 0. _d 0
320 tmpscal2itd(i,j) = 0. _d 0
321 tmpscal3itd(i,j) = 0. _d 0
322 ENDDO
323 ENDDO
324 DO IT=1,SEAICE_multDim
325 DO J=1,sNy
326 DO I=1,sNx
327 C TICES was changed above (item 2), now update TICE as ice volume
328 C weighted average of TICES
329 tmpscal1itd(i,j)=tmpscal1itd(i,j)
330 & + TICES(I,J,IT,bi,bj) * HEFFITD(I,J,IT,bi,bj)
331 tmpscal2itd(i,j)=tmpscal2itd(i,j) + HEFFITD(I,J,IT,bi,bj)
332 C also compute total of AREAITD for diagnostics and SItrArea
333 tmpscal3itd(i,j)=tmpscal3itd(i,j) + AREAITD(I,J,IT,bi,bj)
334 ENDDO
335 ENDDO
336 ENDDO
337 DO J=1,sNy
338 DO I=1,sNx
339 C save pre-ridging ice concentration for diagnostics:
340 C these lines are executed before "ridging" is applied to AREA
341 C hence we execute them here before SEAICE_ITD_REDIST is called
342 C although this means that AREA has not been completely regularized
343 #ifdef ALLOW_DIAGNOSTICS
344 DIAGarrayA(I,J) = tmpscal3itd(i,j)
345 #endif
346 #ifdef ALLOW_SITRACER
347 SItrAREA(I,J,bi,bj,1)=tmpscal3itd(i,j)
348 #endif
349 ENDDO
350 ENDDO
351 C ridge ice according to Lipscomb et al. (2007), Bitz et al. (2001)
352 C Thorndyke et al. (1975), Hibler (1980)
353 CALL SEAICE_DO_RIDGING( bi, bj, myTime, myIter, myThid )
354 C check that all ice thickness categories meet their limits
355 C (includes Hibler-type ridging)
356 CALL SEAICE_ITD_REDIST( bi, bj, myTime, myIter, myThid )
357 C update mean thicknesses HEFF and HSNOW and total ice
358 C concentration AREA to match single category values
359 CALL SEAICE_ITD_SUM ( bi, bj, myTime, myIter, myThid )
360
361 #else /* ifndef SEAICE_ITD */
362
363 #ifdef ALLOW_AUTODIFF_TAMC
364 CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
365 #endif /* ALLOW_AUTODIFF_TAMC */
366 DO J=1,sNy
367 DO I=1,sNx
368 C save pre-ridging ice concentration for diagnostics
369 #ifdef ALLOW_DIAGNOSTICS
370 DIAGarrayA(I,J) = AREA(I,J,bi,bj)
371 #endif /* ALLOW_DIAGNOSTICS */
372 #ifdef ALLOW_SITRACER
373 SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
374 #endif /* ALLOW_SITRACER */
375 C this is the simple Hibler (1979)-type ridging (capping of
376 C concentrations > 1) for the non-ITD sea ice model
377 AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
378 ENDDO
379 ENDDO
380
381 #endif /* SEAICE_ITD */
382
383 #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
384 C end SEAICEadjMODE.EQ.0 statement:
385 ENDIF
386 #endif
387
388 #ifdef ALLOW_DIAGNOSTICS
389 IF ( useDiagnostics ) THEN
390 CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
391 ENDIF
392 #endif /* ALLOW_DIAGNOSTICS */
393
394 C close bi,bj loops
395 ENDDO
396 ENDDO
397
398 #endif /* ALLOW_SEAICE */
399 RETURN
400 END

  ViewVC Help
Powered by ViewVC 1.1.22