/[MITgcm]/MITgcm_contrib/lab_sea_test/growth.F
ViewVC logotype

Contents of /MITgcm_contrib/lab_sea_test/growth.F

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


Revision 1.1 - (show annotations) (download)
Sun Jul 11 06:19:16 2004 UTC (21 years ago) by dimitri
Branch: MAIN
Added growth.F

1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/growth.F,v 1.19 2004/07/06 06:35:04 dimitri Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5
6 CStartOfInterface
7 SUBROUTINE growth( myTime, myIter, myThid )
8 C /==========================================================\
9 C | SUBROUTINE growth |
10 C | o Updata ice thickness and snow depth |
11 C |==========================================================|
12 C \==========================================================/
13 IMPLICIT NONE
14
15 C === Global variables ===
16 #include "SIZE.h"
17 #include "EEPARAMS.h"
18 #include "PARAMS.h"
19 #include "DYNVARS.h"
20 #include "GRID.h"
21 #include "FFIELDS.h"
22 #include "SEAICE_PARAMS.h"
23 #include "SEAICE.h"
24 #include "SEAICE_FFIELDS.h"
25
26 #ifdef ALLOW_AUTODIFF_TAMC
27 # include "tamc.h"
28 #endif
29 C === Routine arguments ===
30 C myTime - Simulation time
31 C myIter - Simulation timestep number
32 C myThid - Thread no. that called this routine.
33 _RL myTime
34 INTEGER myIter, myThid
35 CEndOfInterface
36
37 #ifdef ALLOW_SEAICE
38
39 C === Local variables ===
40 C i,j,bi,bj - Loop counters
41
42 INTEGER i, j, bi, bj
43 _RL TBC, salinity_ice, SDF, Q0, QS
44 _RL GAREA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
45 _RL GHEFF( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
46 _RL AR ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
47
48 C number of surface interface layer
49 INTEGER kSurface
50
51 if ( buoyancyRelation .eq. 'OCEANICP' ) then
52 kSurface = Nr
53 else
54 kSurface = 1
55 endif
56
57 salinity_ice=4.0 _d 0 ! ICE SALINITY
58 TBC=271.2 _d 0-273.16 _d 0 ! FREEZING TEMP. OF SEA WATER
59 SDF=1000.0 _d 0/330.0 _d 0 ! RATIO OF WATER DESITY AND SNOW DENSITY
60 Q0=1.0D-06/302.0 _d +00 ! INVERSE HEAT OF FUSION OF ICE
61 QS=1.1 _d +08 ! HEAT OF FUSION OF SNOW
62
63 DO bj=myByLo(myThid),myByHi(myThid)
64 DO bi=myBxLo(myThid),myBxHi(myThid)
65 c
66 cph(
67 #ifdef ALLOW_AUTODIFF_TAMC
68 act1 = bi - myBxLo(myThid)
69 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
70 act2 = bj - myByLo(myThid)
71 max2 = myByHi(myThid) - myByLo(myThid) + 1
72 act3 = myThid - 1
73 max3 = nTx*nTy
74 act4 = ikey_dynamics - 1
75 iicekey = (act1 + 1) + act2*max1
76 & + act3*max1*max2
77 & + act4*max1*max2*max3
78 #endif /* ALLOW_AUTODIFF_TAMC */
79 c
80 #ifdef ALLOW_AUTODIFF_TAMC
81 CADJ STORE theta(:,:,:,bi,bj)= comlev1_bibj,
82 CADJ & key = iicekey, byte = isbyte
83 CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
84 CADJ & key = iicekey, byte = isbyte
85 CADJ STORE atemp(:,:,bi,bj) = comlev1_bibj,
86 CADJ & key = iicekey, byte = isbyte
87 #endif /* ALLOW_AUTODIFF_TAMC */
88 cph)
89 DO J=1,sNy
90 DO I=1,sNx
91 SEAICE_SALT(I,J,bi,bj)=ZERO
92 ENDDO
93 ENDDO
94 #ifdef ALLOW_AUTODIFF_TAMC
95 CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
96 CADJ & key = iicekey, byte = isbyte
97 CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
98 CADJ & key = iicekey, byte = isbyte
99 #endif /* ALLOW_AUTODIFF_TAMC */
100 DO J=1,sNy
101 DO I=1,sNx
102 AR(I,J,bi,bj)=MIN(AREA(I,J,2,bi,bj),
103 & HEFF(I,J,2,bi,bj)*1.0 _d +04)
104 ENDDO
105 ENDDO
106 #ifdef ALLOW_AUTODIFF_TAMC
107 CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
108 CADJ & key = iicekey, byte = isbyte
109 #endif /* ALLOW_AUTODIFF_TAMC */
110 DO J=1,sNy
111 DO I=1,sNx
112 C-- Create or melt sea-ice so that first-level oceanic temperature
113 C is approximately at the freezing point when there is sea-ice.
114 C Initially the units of YNEG are m of sea-ice.
115 C The factor dRf(1)/72.0764, used to convert temperature
116 C change in deg K to m of sea-ice, is approximately:
117 C dRf(1) * (sea water heat capacity = 3996 J/kg/K)
118 C * (density of sea-water = 1026 kg/m^3)
119 C / (latent heat of fusion of sea-ice = 334000 J/kg)
120 C / (density of sea-ice = 910 kg/m^3)
121 C Negative YNEG leads to ice growth.
122 C Positive YNEG leads to ice melting.
123 YNEG(I,J,bi,bj)=(theta(I,J,1,bi,bj)-TBC)*.01
124 & *dRf(1)/72.0764 _d 0
125 GHEFF(I,J)=HEFF(I,J,1,bi,bj)
126 HEFF(I,J,1,bi,bj)=MAX(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj))
127 YNEG(I,J,bi,bj)=GHEFF(I,J)-HEFF(I,J,1,bi,bj)
128 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)-YNEG(I,J,bi,bj)
129 C-- Now convert YNEG back to deg K.
130 YNEG(I,J,bi,bj)=YNEG(I,J,bi,bj)*recip_dRf(1)*72.0764 _d 0
131 ENDDO
132 ENDDO
133 c
134 ENDDO
135 ENDDO
136
137 cph(
138 #ifdef ALLOW_AUTODIFF_TAMC
139 CADJ STORE area = comlev1, key = ikey_dynamics
140 CADJ STORE atemp = comlev1, key = ikey_dynamics
141 CADJ STORE heff = comlev1, key = ikey_dynamics
142 CADJ STORE hsnow = comlev1, key = ikey_dynamics
143 CADJ STORE lwdown = comlev1, key = ikey_dynamics
144 CADJ STORE tice = comlev1, key = ikey_dynamics
145 CADJ STORE uwind = comlev1, key = ikey_dynamics
146 CADJ STORE vwind = comlev1, key = ikey_dynamics
147 # ifdef SEAICE_MULTILEVEL
148 CADJ STORE tices = comlev1, key = ikey_dynamics
149 # endif
150 #endif /* ALLOW_AUTODIFF_TAMC */
151 cph)
152 C GROWTH SUBROUTINE CALCULATES TOTAL GROWTH TENDENCIES,
153 C INCLUDING SNOWFALL
154 CALL GROATB(A22,myThid)
155
156 DO bj=myByLo(myThid),myByHi(myThid)
157 DO bi=myBxLo(myThid),myBxHi(myThid)
158 cph(
159 #ifdef ALLOW_AUTODIFF_TAMC
160 act1 = bi - myBxLo(myThid)
161 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
162 act2 = bj - myByLo(myThid)
163 max2 = myByHi(myThid) - myByLo(myThid) + 1
164 act3 = myThid - 1
165 max3 = nTx*nTy
166 act4 = ikey_dynamics - 1
167 iicekey = (act1 + 1) + act2*max1
168 & + act3*max1*max2
169 & + act4*max1*max2*max3
170 #endif /* ALLOW_AUTODIFF_TAMC */
171 c
172 #ifdef ALLOW_AUTODIFF_TAMC
173 CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
174 CADJ & key = iicekey, byte = isbyte
175 CADJ STORE heff(:,:,:,bi,bj) = comlev1_bibj,
176 CADJ & key = iicekey, byte = isbyte
177 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
178 CADJ & key = iicekey, byte = isbyte
179 CADJ STORE fo(:,:,bi,bj) = comlev1_bibj,
180 CADJ & key = iicekey, byte = isbyte
181 CADJ STORE fice(:,:,bi,bj) = comlev1_bibj,
182 CADJ & key = iicekey, byte = isbyte
183 #endif /* ALLOW_AUTODIFF_TAMC */
184 cph)
185 C NOW CALCULATE CORRECTED GROWTH
186 DO J=1,sNy
187 DO I=1,sNx
188 GHEFF(I,J)=-SEAICE_deltaTtherm*FICE(I,J,bi,bj)
189 GAREA(I,J)=HSNOW(I,J,bi,bj)*QS
190 IF(GHEFF(I,J).GT.ZERO.AND.GHEFF(I,J).LE.GAREA(I,J)) THEN
191 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS
192 C SNOW CONVERTED INTO WATER AND THEN INTO ICE
193 C The factor 0.920 is used to convert m of sea-ice to m of freshwater
194 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
195 & -(GHEFF(I,J)/QS)/SDF/0.920 _d 0*AR(I,J,bi,bj)
196 FICE(I,J,bi,bj)=ZERO
197 ELSE IF(GHEFF(I,J).GT.GAREA(I,J)) THEN
198 FICE(I,J,bi,bj)=-(GHEFF(I,J)-GAREA(I,J))/SEAICE_deltaTtherm
199 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
200 & -HSNOW(I,J,bi,bj)/SDF/0.920 _d 0*AR(I,J,bi,bj)
201 HSNOW(I,J,bi,bj)=0.0
202 END IF
203
204 ENDDO
205 ENDDO
206
207 C NOW GET TOTAL GROWTH RATE
208 DO J=1,sNy
209 DO I=1,sNx
210 FHEFF(I,J,bi,bj)=FICE(I,J,bi,bj)*AR(I,J,bi,bj)
211 & +(ONE-AR(I,J,bi,bj))*FO(I,J,bi,bj)
212 ENDDO
213 ENDDO
214
215
216 C NOW UPDATE AREA
217 DO J=1,sNy
218 DO I=1,sNx
219 GHEFF(I,J)=-SEAICE_deltaTtherm*FHEFF(I,J,bi,bj)*Q0
220 GAREA(I,J)=SEAICE_deltaTtherm*FO(I,J,bi,bj)*Q0
221 GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))
222 GAREA(I,J)=MAX(ZERO,GAREA(I,J))
223 HCORR(I,J,bi,bj)=MIN(ZERO,GHEFF(I,J))
224 ENDDO
225 ENDDO
226 DO J=1,sNy
227 DO I=1,sNx
228 GAREA(I,J)=TWO*(ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO
229 & +HALF*HCORR(I,J,bi,bj)*AREA(I,J,2,bi,bj)
230 & /(HEFF(I,J,1,bi,bj)+.00001 _d 0)
231 AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+GAREA(I,J)
232 ENDDO
233 ENDDO
234
235 C NOW UPDATE HEFF
236 DO J=1,sNy
237 DO I=1,sNx
238 GHEFF(I,J)=-SEAICE_deltaTtherm*
239 & FICE(I,J,bi,bj)*Q0*AR(I,J,bi,bj)
240 GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))
241 HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)+GHEFF(I,J)
242 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)+GHEFF(I,J)
243 C NOW CALCULATE QNETI UNDER ICE IF ANY
244 QNETI(I,J,bi,bj)=(GHEFF(I,J)-SEAICE_deltaTtherm*
245 & FICE(I,J,bi,bj)*Q0*AR(I,J,bi,bj))/Q0/SEAICE_deltaTtherm
246 ENDDO
247 ENDDO
248
249 C NOW GET TOTAL QNET AND QSW
250 DO J=1,sNy
251 DO I=1,sNx
252 QNET(I,J,bi,bj)=QNETI(I,J,bi,bj)*AR(I,J,bi,bj)
253 & +(ONE-AR(I,J,bi,bj))*QNETO(I,J,bi,bj)
254 QSW(I,J,bi,bj)=QSWI(I,J,bi,bj)*AR(I,J,bi,bj)
255 & +(ONE-AR(I,J,bi,bj))*QSWO(I,J,bi,bj)
256 #ifndef SHORTWAVE_HEATING
257 QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+QSW(I,J,bi,bj)
258 #endif
259 C Add YNEG contribution to QNET
260 QNET(I,J,bi,bj)=QNET(I,J,bi,bj)
261 & +YNEG(I,J,bi,bj)/SEAICE_deltaTtherm*maskC(I,J,1,bi,bj)
262 & *HeatCapacity_Cp*recip_horiVertRatio*rhoConst
263 & *drF(kSurface)*hFacC(i,j,kSurface,bi,bj)
264 ENDDO
265 ENDDO
266
267 C NOW UPDATE OTHER THINGS
268 DO J=1,sNy
269 DO I=1,sNx
270 IF(FICE(I,J,bi,bj).GT.ZERO) THEN
271 C FREEZING, PRECIP ADDED AS SNOW
272 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+SEAICE_deltaTtherm*
273 & PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
274 ELSE
275 C ADD PRECIP AS RAIN, WATER CONVERTED INTO ICE BY /0.920 _d 0
276 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
277 & -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*
278 & SEAICE_deltaTtherm/0.920 _d 0
279 ENDIF
280 c Now add in precip over open water directly into ocean as negative salt
281 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
282 & -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
283 & *SEAICE_deltaTtherm/0.920 _d 0
284 C NOW GET FRESH WATER FLUX
285 EmPmR(I,J,bi,bj)= maskC(I,J,1,bi,bj)*(
286 & EVAP(I,J,bi,bj)-RUNOFF(I,J,bi,bj)
287 & +SEAICE_SALT(I,J,bi,bj)*0.92 _d 0/SEAICE_deltaTtherm
288 & )
289 ENDDO
290 ENDDO
291
292 #ifdef SEAICE_DEBUG
293 c CALL PLOT_FIELD_XYRS( UWIND,'Current UWIND ', myIter, myThid )
294 c CALL PLOT_FIELD_XYRS( VWIND,'Current VWIND ', myIter, myThid )
295 CALL PLOT_FIELD_XYRS( GWATX,'Current GWATX ', myIter, myThid )
296 CALL PLOT_FIELD_XYRS( GWATY,'Current GWATY ', myIter, myThid )
297 CALL PLOT_FIELD_XYRL( FO,'Current FO ', myIter, myThid )
298 CALL PLOT_FIELD_XYRL( FHEFF,'Current FHEFF ', myIter, myThid )
299 CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
300 CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
301 CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
302 DO j=1-OLy,sNy+OLy
303 DO i=1-OLx,sNx+OLx
304 GHEFF(I,J)=SQRT(UICE(I,J,1,bi,bj)**2+VICE(I,J,1,bi,bj)**2)
305 GAREA(I,J)=HEFF(I,J,1,bi,bj)
306 print*,'I J QNET:',I, J, QNET(i,j,bi,bj), QSW(I,J,bi,bj)
307 ENDDO
308 ENDDO
309 CALL PLOT_FIELD_XYRL( GHEFF,'Current UICE ', myIter, myThid )
310 CALL PLOT_FIELD_XYRL( GAREA,'Current HEFF ', myIter, myThid )
311 DO j=1-OLy,sNy+OLy
312 DO i=1-OLx,sNx+OLx
313 if(HEFF(i,j,1,bi,bj).gt.1.) then
314 print '(A,2i4,3f10.2)','#### i j heff theta yneg',i,j,
315 & HEFF(i,j,1,bi,bj),theta(I,J,1,bi,bj),yneg(I,J,bi,bj)
316 print '(A,3f10.2)','QSW, QNET before/after correction',
317 & QSW(I,J,bi,bj),QNETI(I,J,bi,bj)*AR(I,J,bi,bj)
318 & +(ONE-AR(I,J,bi,bj))*QNETO(I,J,bi,bj), QNET(I,J,bi,bj)
319 endif
320 ENDDO
321 ENDDO
322 #endif /* SEAICE_DEBUG */
323
324 crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ?
325 #define DO_WE_NEED_THIS
326 C NOW ZERO OUTSIDE POINTS
327 DO J=1,sNy
328 DO I=1,sNx
329 C NOW SET AREA(I,J,1,bi,bj)=0 WHERE NO ICE IS
330 AREA(I,J,1,bi,bj)=MIN(AREA(I,J,1,bi,bj)
331 & ,HEFF(I,J,1,bi,bj)/.0001 _d 0)
332 ENDDO
333 ENDDO
334 #ifdef ALLOW_AUTODIFF_TAMC
335 CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
336 CADJ & key = iicekey, byte = isbyte
337 #endif /* ALLOW_AUTODIFF_TAMC */
338 DO J=1,sNy
339 DO I=1,sNx
340 C NOW TRUNCATE AREA
341 #ifdef DO_WE_NEED_THIS
342 AREA(I,J,1,bi,bj)=MIN(ONE,AREA(I,J,1,bi,bj))
343 ENDDO
344 ENDDO
345 #ifdef ALLOW_AUTODIFF_TAMC
346 CADJ STORE area(:,:,:,bi,bj) = comlev1_bibj,
347 CADJ & key = iicekey, byte = isbyte
348 #endif /* ALLOW_AUTODIFF_TAMC */
349 DO J=1,sNy
350 DO I=1,sNx
351 AREA(I,J,1,bi,bj)=MAX(ZERO,AREA(I,J,1,bi,bj))
352 HSNOW(I,J,bi,bj)=MAX(ZERO,HSNOW(I,J,bi,bj))
353 #endif
354 AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
355 HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
356 #ifdef DO_WE_NEED_THIS
357 c HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))
358 #endif
359 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
360 ENDDO
361 ENDDO
362
363 ENDDO
364 ENDDO
365
366 #endif /* ALLOW_SEAICE */
367
368 RETURN
369 END

  ViewVC Help
Powered by ViewVC 1.1.22