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

Annotation of /MITgcm_contrib/lab_sea_test/growth.F

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


Revision 1.2 - (hide annotations) (download)
Mon Jul 12 01:00:20 2004 UTC (19 years, 10 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 dimitri 1.2 C $Header: /u/gcmpack/MITgcm_contrib/lab_sea_test/growth.F,v 1.1 2004/07/11 06:19:16 dimitri Exp $
2 dimitri 1.1 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 dimitri 1.2 _RL mymin_R8, mymax_R8
52     external mymin_R8, mymax_R8
53    
54 dimitri 1.1 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 dimitri 1.2 AR(I,J,bi,bj)=MYMIN_R8(AREA(I,J,2,bi,bj),
106 dimitri 1.1 & 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 dimitri 1.2 HEFF(I,J,1,bi,bj)=MYMAX_R8(ZERO,HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj))
130 dimitri 1.1 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 dimitri 1.2 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 dimitri 1.1 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 dimitri 1.2 GHEFF(I,J)=-ONE*MYMIN_R8(HEFF(I,J,1,bi,bj),GHEFF(I,J))
244 dimitri 1.1 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 dimitri 1.2 AREA(I,J,1,bi,bj)=MYMIN_R8(AREA(I,J,1,bi,bj)
334 dimitri 1.1 & ,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 dimitri 1.2 AREA(I,J,1,bi,bj)=MYMIN_R8(ONE,AREA(I,J,1,bi,bj))
346 dimitri 1.1 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 dimitri 1.2 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 dimitri 1.1 #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 dimitri 1.2 c HEFF(I,J,1,bi,bj)=MYMIN_R8(MAX_HEFF,HEFF(I,J,1,bi,bj))
361 dimitri 1.1 #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