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 |