/[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.2 - (show annotations) (download)
Mon Jul 12 01:00:20 2004 UTC (19 years, 9 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +15 -12 lines
added my_min_max for pkg/seaice routines

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

  ViewVC Help
Powered by ViewVC 1.1.22