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

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

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

Revision 1.5 - (show annotations) (download)
Tue Feb 18 05:33:55 2003 UTC (21 years, 2 months ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint48f_post, checkpoint48i_post, checkpoint50e_post, checkpoint50c_post, checkpoint48h_post, checkpoint50c_pre, checkpoint50d_pre, checkpoint50b_pre, checkpoint49, checkpoint48g_post, checkpoint50, checkpoint50d_post, checkpoint50b_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint50e_pre
Changes since 1.4: +15 -6 lines 
Merging from release1_p12:
o Modifications for using pkg/exf with pkg/seaice
  - improved description of the various forcing configurations
  - added basic radiation bulk formulae to pkg/exf
  - units/sign fix for evap computation in exf_getffields.F
  - updated verification/global_with_exf/results/output.txt
o Added pkg/sbo for computing IERS Special Bureau for the Oceans
  (SBO) core products, including oceanic mass, center-of-mass,
  angular, and bottom pressure (see pkg/sbo/README.sbo).
o Lower bound for viscosity/diffusivity in pkg/kpp/kpp_routines.F
  to avoid negative values in shallow regions.
  - updated verification/natl_box/results/output.txt
  - updated verification/lab_sea/results/output.txt
o MPI gather, scatter: eesupp/src/gather_2d.F and scatter_2d.F
o Added useSingleCpuIO option (see PARAMS.h).
o Updated useSingleCpuIO option in mdsio_writefield.F to
  work with multi-field files, e.g., for single-file pickup.
o pkg/seaice:
  - bug fix in growth.F: QNET for no shortwave case
  - added HeffFile for specifying initial sea-ice thickness
  - changed SEAICE_EXTERNAL_FLUXES wind stress implementation
o Added missing /* */ to CPP comments in pkg/seaice, pkg/exf,
  kpp_transport_t.F, forward_step.F, and the_main_loop.F
o pkg/seaice:
  - adjoint-friendly modifications
  - added a SEAICE_WRITE_PICKUP at end of the_model_main.F
1 C $Header:
2
3 #include "SEAICE_OPTIONS.h"
4
5 CStartOfInterface
6 SUBROUTINE growth( myTime, myIter, myThid )
7 C /==========================================================\
8 C | SUBROUTINE growth |
9 C | o Updata ice thickness and snow depth |
10 C |==========================================================|
11 C \==========================================================/
12 IMPLICIT NONE
13
14 C === Global variables ===
15 #include "SIZE.h"
16 #include "EEPARAMS.h"
17 #include "PARAMS.h"
18 #include "DYNVARS.h"
19 #include "GRID.h"
20 #include "FFIELDS.h"
21 #include "SEAICE_PARAMS.h"
22 #include "SEAICE.h"
23 #include "SEAICE_FFIELDS.h"
24 #include "SEAICE_EXTERNAL.h"
25
26 C === Routine arguments ===
27 C myTime - Simulation time
28 C myIter - Simulation timestep number
29 C myThid - Thread no. that called this routine.
30 _RL myTime
31 INTEGER myIter, myThid
32 CEndOfInterface
33
34 #ifdef ALLOW_SEAICE
35
36 C === Local variables ===
37 C i,j,bi,bj - Loop counters
38
39 INTEGER i, j, bi, bj
40 _RL TBC, salinity_ice, SDF, Q0, QS
41 _RL theta_old
42
43 _RL GAREA( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
44 _RL GHEFF( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy )
45 _RL AR ( 1-OLx:sNx+OLx, 1-OLy:sNy+OLy, nSx, nSy )
46
47 salinity_ice=4.0 _d 0 ! ICE SALINITY
48 TBC=271.2 _d 0-273.16 _d 0 ! FREEZING TEMP. OF SEA WATER
49 SDF=1000.0 _d 0/330.0 _d 0 ! RATIO OF WATER DESITY AND SNOW DENSITY
50 Q0=1.0D-06/302.0 _d +00 ! INVERSE HEAT OF FUSION OF ICE
51 QS=1.1 _d +08 ! HEAT OF FUSION OF SNOW
52
53 DO bj=myByLo(myThid),myByHi(myThid)
54 DO bi=myBxLo(myThid),myBxHi(myThid)
55 DO J=1,sNy
56 DO I=1,sNx
57 SEAICE_SALT(I,J,bi,bj)=ZERO
58 WATR(I,J,bi,bj)=ZERO
59 AR(I,J,bi,bj)=MIN(AREA(I,J,2,bi,bj),
60 & HEFF(I,J,2,bi,bj)*1.0 _d +04)
61 C NOW BALANCE THE HEAT IN OCEAN FIRT LEVEL
62 C Here the units of YNEG are m of sea-ice.
63 C The factor dRf(1)/72.0764, used to convert temperature
64 C change in deg K to m of sea-ice, is approximately:
65 C dRf(1) * (sea water heat capacity = 3996 J/kg/K)
66 C * (density of sea-water = 1026 kg/m^3)
67 C / (latent heat of fusion of sea-ice = 334000 J/kg)
68 C / (density of sea-ice = 910 kg/m^3)
69 YNEG(I,J,bi,bj)=(theta(I,J,1,bi,bj)-TBC)*dRf(1)/72.0764 _d 0
70 IF(YNEG(I,J,bi,bj).LE.ZERO) THEN
71 C SUPERCOOLING, CONVERT TO ICE
72 HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj)
73 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)-YNEG(I,J,bi,bj)
74 C HERE YNEG IS THE OCEAN TEMP. DIFFERENCE WITH SUPER-COOLING
75 YNEG(I,J,bi,bj)=TBC-theta(I,J,1,bi,bj)
76 theta(I,J,1,bi,bj)=TBC
77 ELSE
78 C NOW CORRECT ICE THICKNESS
79 GHEFF(I,J)=HEFF(I,J,1,bi,bj)
80 HEFF(I,J,3,bi,bj)=HEFF(I,J,1,bi,bj)-YNEG(I,J,bi,bj)
81 HEFF(I,J,1,bi,bj)=MAX(ZERO,HEFF(I,J,3,bi,bj))
82 C NOW CORRECT YNEG
83 YNEG(I,J,bi,bj)=HEFF(I,J,1,bi,bj)-HEFF(I,J,3,bi,bj)
84 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
85 & +(HEFF(I,J,1,bi,bj)-GHEFF(I,J))
86 theta_old=theta(I,J,1,bi,bj)
87 theta(I,J,1,bi,bj)=
88 & (YNEG(I,J,bi,bj)*72.0764 _d 0/dRf(1)+TBC)
89 C HERE YNEG IS THE OCEAN TEMP. DIFFERENCE WITH NO SUPER-COOLING
90 YNEG(I,J,bi,bj)=theta(I,J,1,bi,bj)-theta_old
91 ENDIF
92 C surfaceTendencyTice is the equivalent amount of surface heating
93 C implied by above change in surface level temperature (YNEG).
94 C Units are degrees/s (>0 for ocean warming).
95 surfaceTendencyTice(I,J,bi,bj)=YNEG(I,J,bi,bj)/deltaTtracer
96 ENDDO
97 ENDDO
98 ENDDO
99 ENDDO
100
101 C GROWTH SUBROUTINE CALCULATES TOTAL GROWTH TENDENCIES,
102 C INCLUDING SNOWFALL
103 CALL GROATB(A22,myThid)
104
105 DO bj=myByLo(myThid),myByHi(myThid)
106 DO bi=myBxLo(myThid),myBxHi(myThid)
107 C NOW CALCULATE CORRECTED GROWTH
108 DO J=1,sNy
109 DO I=1,sNx
110 GHEFF(I,J)=-DELTAT*FICE(I,J,bi,bj)
111 GAREA(I,J)=HSNOW(I,J,bi,bj)*QS
112 IF(GHEFF(I,J).GT.ZERO.AND.GHEFF(I,J).LE.GAREA(I,J)) THEN
113 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)-GHEFF(I,J)/QS
114 C SNOW CONVERTED INTO WATER AND THEN INTO ICE
115 C The factor 0.920 is used to convert m of sea-ice to m of freshwater
116 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
117 & -(GHEFF(I,J)/QS)/SDF/0.920 _d 0*AR(I,J,bi,bj)
118 WATR(I,J,bi,bj)=WATR(I,J,bi,bj)
119 & -(GHEFF(I,J)/QS)/SDF/0.920 _d 0*AR(I,J,bi,bj)
120 FICE(I,J,bi,bj)=ZERO
121 ELSE IF(GHEFF(I,J).GT.GAREA(I,J)) THEN
122 FICE(I,J,bi,bj)=-(GHEFF(I,J)-GAREA(I,J))/DELTAT
123 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
124 & -HSNOW(I,J,bi,bj)/SDF/0.920 _d 0*AR(I,J,bi,bj)
125 WATR(I,J,bi,bj)=WATR(I,J,bi,bj)
126 & -HSNOW(I,J,bi,bj)/SDF/0.920*AR(I,J,bi,bj)
127 HSNOW(I,J,bi,bj)=0.0
128 END IF
129
130 ENDDO
131 ENDDO
132
133 C NOW GET TOTAL GROWTH RATE
134 DO J=1,sNy
135 DO I=1,sNx
136 FHEFF(I,J,bi,bj)=FICE(I,J,bi,bj)*AR(I,J,bi,bj)
137 & +(ONE-AR(I,J,bi,bj))*FO(I,J,bi,bj)
138 ENDDO
139 ENDDO
140
141
142 C NOW UPDATE AREA
143 DO J=1,sNy
144 DO I=1,sNx
145 GHEFF(I,J)=-DELTAT*FHEFF(I,J,bi,bj)*Q0
146 GAREA(I,J)=DELTAT*FO(I,J,bi,bj)*Q0
147 GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))
148 GAREA(I,J)=MAX(ZERO,GAREA(I,J))
149 HCORR(I,J,bi,bj)=MIN(ZERO,GHEFF(I,J))
150 ENDDO
151 ENDDO
152 DO J=1,sNy
153 DO I=1,sNx
154 GAREA(I,J)=TWO*(ONE-AREA(I,J,2,bi,bj))*GAREA(I,J)/HO
155 & +HALF*HCORR(I,J,bi,bj)*AREA(I,J,2,bi,bj)
156 & /(HEFF(I,J,1,bi,bj)+.00001 _d 0)
157 AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)+GAREA(I,J)
158 ENDDO
159 ENDDO
160
161 C NOW UPDATE HEFF
162 DO J=1,sNy
163 DO I=1,sNx
164 GHEFF(I,J)=-DELTAT*FICE(I,J,bi,bj)*Q0*AR(I,J,bi,bj)
165 GHEFF(I,J)=-ONE*MIN(HEFF(I,J,1,bi,bj),GHEFF(I,J))
166 HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)+GHEFF(I,J)
167 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)+GHEFF(I,J)
168 C NOW CALCULATE QNETI UNDER ICE IF ANY
169 QNETI(I,J,bi,bj)=(GHEFF(I,J)-DELTAT*FICE(I,J,bi,bj)
170 & *Q0*AR(I,J,bi,bj))/Q0/DELTAT
171 ENDDO
172 ENDDO
173
174 C NOW GET TOTAL QNET AND QSW
175 DO J=1,sNy
176 DO I=1,sNx
177 QNET(I,J,bi,bj)=QNETI(I,J,bi,bj)*AR(I,J,bi,bj)
178 & +(ONE-AR(I,J,bi,bj))*QNETO(I,J,bi,bj)
179 QSW(I,J,bi,bj)=QSWI(I,J,bi,bj)*AR(I,J,bi,bj)
180 & +(ONE-AR(I,J,bi,bj))*QSWO(I,J,bi,bj)
181 #ifndef SHORTWAVE_HEATING
182 QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+QSW(I,J,bi,bj)
183 #endif
184 ENDDO
185 ENDDO
186
187 C NOW UPDATE OTHER THINGS
188 DO J=1,sNy
189 DO I=1,sNx
190 IF(FICE(I,J,bi,bj).GT.ZERO) THEN
191 C FREEZING, PRECIP ADDED AS SNOW
192 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)
193 & +DELTAT*PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*SDF
194 ELSE
195 C ADD PRECIP AS RAIN, WATER CONVERTED INTO ICE BY /0.920 _d 0
196 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
197 & -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*DELTAT/0.920 _d 0
198 WATR(I,J,bi,bj)=WATR(I,J,bi,bj)
199 & -PRECIP(I,J,bi,bj)*AREA(I,J,2,bi,bj)*DELTAT/0.920 _d 0
200 ENDIF
201 c Now add in precip over open water directly into ocean as negative salt
202 SEAICE_SALT(I,J,bi,bj)=SEAICE_SALT(I,J,bi,bj)
203 & -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
204 & *DELTAT/0.920 _d 0
205 WATR(I,J,bi,bj)=WATR(I,J,bi,bj)
206 & -PRECIP(I,J,bi,bj)*(ONE-AREA(I,J,2,bi,bj))
207 & *DELTAT/0.920 _d 0
208 C NOW GET FRESH WATER FLUX
209 EmPmR(I,J,bi,bj)=EVAP(I,J,bi,bj)-RUNOFF(I,J,bi,bj)
210 & +SEAICE_SALT(I,J,bi,bj)*0.92 _d 0/DELTAT
211 ENDDO
212 ENDDO
213
214 #ifdef SEAICE_DEBUG
215 c CALL PLOT_FIELD_XYRS( UWIND,'Current UWIND ', myIter, myThid )
216 c CALL PLOT_FIELD_XYRS( VWIND,'Current VWIND ', myIter, myThid )
217 CALL PLOT_FIELD_XYRS( GWATX,'Current GWATX ', myIter, myThid )
218 CALL PLOT_FIELD_XYRS( GWATY,'Current GWATY ', myIter, myThid )
219 CALL PLOT_FIELD_XYRL( FO,'Current FO ', myIter, myThid )
220 CALL PLOT_FIELD_XYRL( FHEFF,'Current FHEFF ', myIter, myThid )
221 CALL PLOT_FIELD_XYRL( QSW,'Current QSW ', myIter, myThid )
222 CALL PLOT_FIELD_XYRL( QNET,'Current QNET ', myIter, myThid )
223 CALL PLOT_FIELD_XYRL( EmPmR,'Current EmPmR ', myIter, myThid )
224 DO j=1-OLy,sNy+OLy
225 DO i=1-OLx,sNx+OLx
226 GHEFF(I,J)=SQRT(UICE(I,J,1,bi,bj)**2+VICE(I,J,1,bi,bj)**2)
227 GAREA(I,J)=HEFF(I,J,1,bi,bj)
228 print*,'I J QNET:',I, J, QNET(i,j,bi,bj), QSW(I,J,bi,bj)
229 ENDDO
230 ENDDO
231 CALL PLOT_FIELD_XYRL( GHEFF,'Current UICE ', myIter, myThid )
232 CALL PLOT_FIELD_XYRL( GAREA,'Current HEFF ', myIter, myThid )
233 DO j=1-OLy,sNy+OLy
234 DO i=1-OLx,sNx+OLx
235 if(HEFF(i,j,1,bi,bj).gt.1.) then
236 print '(A,2i4,3f10.2)','#### i j heff theta yneg',i,j,
237 & HEFF(i,j,1,bi,bj),theta(I,J,1,bi,bj),yneg(I,J,bi,bj)
238 print '(A,3f10.2)','QSW, QNET before/after correction',
239 & QSW(I,J,bi,bj),QNETI(I,J,bi,bj)*AR(I,J,bi,bj)
240 & +(ONE-AR(I,J,bi,bj))*QNETO(I,J,bi,bj), QNET(I,J,bi,bj)
241 endif
242 ENDDO
243 ENDDO
244 #endif /* SEAICE_DEBUG */
245
246 crg Added by Ralf Giering: do we need DO_WE_NEED_THIS ?
247 #define DO_WE_NEED_THIS
248 C NOW ZERO OUTSIDE POINTS
249 DO J=1,sNy
250 DO I=1,sNx
251 C NOW SET AREA(I,J,1,bi,bj)=0 WHERE NO ICE IS
252 AREA(I,J,1,bi,bj)=MIN(AREA(I,J,1,bi,bj)
253 & ,HEFF(I,J,1,bi,bj)/.0001 _d 0)
254 C NOW TRUNCATE AREA
255 #ifdef DO_WE_NEED_THIS
256 AREA(I,J,1,bi,bj)=MIN(ONE,AREA(I,J,1,bi,bj))
257 AREA(I,J,1,bi,bj)=MAX(ZERO,AREA(I,J,1,bi,bj))
258 HSNOW(I,J,bi,bj)=MAX(ZERO,HSNOW(I,J,bi,bj))
259 #endif
260 AREA(I,J,1,bi,bj)=AREA(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
261 HEFF(I,J,1,bi,bj)=HEFF(I,J,1,bi,bj)*HEFFM(I,J,bi,bj)
262 #ifdef DO_WE_NEED_THIS
263 HEFF(I,J,1,bi,bj)=MIN(MAX_HEFF,HEFF(I,J,1,bi,bj))
264 #endif
265 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
266 ENDDO
267 ENDDO
268
269 ENDDO
270 ENDDO
271
272 #endif /* ALLOW_SEAICE */
273
274 RETURN
275 END
  ViewVC Help
Powered by ViewVC 1.1.22