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

Annotation of /MITgcm/pkg/seaice/seaice_solve4temp.F

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

Revision 1.15 - (hide annotations) (download)
Sun Jun 19 23:01:54 2011 UTC (12 years, 11 months ago) by ifenty
Branch: MAIN
CVS Tags: checkpoint62z
Changes since 1.14: +4 -3 lines 
Removed redundant minimum ice thickness check (hice_tmp = max(HICE_ACTUAL, 0.05)) in
the evolution branch now that seaice_growth does the regularization (SEAICE_hice_reg).
Removing that line doesn't change any verification results.
1 ifenty 1.15 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_solve4temp.F,v 1.14 2011/06/19 02:31:40 ifenty Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: SEAICE_SOLVE4TEMP
8     C !INTERFACE:
9     SUBROUTINE SEAICE_SOLVE4TEMP(
10     I UG, HICE_ACTUAL, HSNOW_ACTUAL,
11     U TSURF,
12     O F_ia, IcePenetSWFlux,
13 mlosch 1.10 O FWsublim,
14 jmc 1.1 I bi, bj, myTime, myIter, myThid )
15    
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | SUBROUTINE SOLVE4TEMP
19     C | o Calculate ice growth rate, surface fluxes and
20     C | temperature of ice surface.
21     C | see Hibler, MWR, 108, 1943-1973, 1980
22     C *==========================================================*
23     C \ev
24    
25     C !USES:
26     IMPLICIT NONE
27     C === Global variables ===
28     #include "SIZE.h"
29     #include "GRID.h"
30     #include "EEPARAMS.h"
31 jmc 1.3 #include "PARAMS.h"
32 jmc 1.1 #include "FFIELDS.h"
33 heimbach 1.13 #include "SEAICE_SIZE.h"
34     #include "SEAICE_PARAMS.h"
35 jmc 1.1 #include "SEAICE.h"
36     #ifdef SEAICE_VARIABLE_FREEZING_POINT
37     #include "DYNVARS.h"
38     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
39     #ifdef ALLOW_EXF
40     # include "EXF_OPTIONS.h"
41     # include "EXF_FIELDS.h"
42     #endif
43 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
44     # include "tamc.h"
45     #endif
46 jmc 1.1
47     C !INPUT/OUTPUT PARAMETERS
48     C === Routine arguments ===
49     C INPUT:
50     C UG :: thermal wind of atmosphere
51     C HICE_ACTUAL :: actual ice thickness
52     C HSNOW_ACTUAL :: actual snow thickness
53     C TSURF :: surface temperature of ice in Kelvin, updated
54     C bi,bj :: loop indices
55     C OUTPUT:
56 jmc 1.12 C F_io_net :: net upward conductive heat flux through ice at the base
57 mlosch 1.10 C of the ice
58 jmc 1.1 C F_ia_net :: net heat flux divergence at the sea ice/snow surface:
59     C includes ice conductive fluxes and atmospheric fluxes (W/m^2)
60     C F_ia :: upward sea ice/snow surface heat flux to atmosphere (W/m^2)
61     C IcePenetSWFlux :: short wave heat flux under ice
62 jmc 1.12 C FWsublim :: fresh water (mass) flux implied by latent heat of
63 mlosch 1.10 C sublimation (kg/m^2/s)
64 jmc 1.1 _RL UG (1:sNx,1:sNy)
65     _RL HICE_ACTUAL (1:sNx,1:sNy)
66     _RL HSNOW_ACTUAL (1:sNx,1:sNy)
67     _RL TSURF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
68 jmc 1.12 c _RL F_io_net (1:sNx,1:sNy)
69     c _RL F_ia_net (1:sNx,1:sNy)
70 jmc 1.1 _RL F_ia (1:sNx,1:sNy)
71     _RL IcePenetSWFlux (1:sNx,1:sNy)
72 mlosch 1.10 _RL FWsublim (1:sNx,1:sNy)
73 jmc 1.1 INTEGER bi, bj
74     _RL myTime
75     INTEGER myIter, myThid
76    
77     C !LOCAL VARIABLES:
78     C === Local variables ===
79 jmc 1.12 _RL F_io_net (1:sNx,1:sNy)
80     _RL F_ia_net (1:sNx,1:sNy)
81 jmc 1.6 #ifndef SEAICE_SOLVE4TEMP_LEGACY
82 jmc 1.1 _RL F_swi (1:sNx,1:sNy)
83     _RL F_lwd (1:sNx,1:sNy)
84     _RL F_lwu (1:sNx,1:sNy)
85     _RL F_sens (1:sNx,1:sNy)
86 gforget 1.11 _RL hice_tmp
87 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
88 mlosch 1.10 _RL F_lh (1:sNx,1:sNy)
89 jmc 1.1 _RL F_c (1:sNx,1:sNy)
90     _RL qhice (1:sNx,1:sNy)
91    
92     _RL AbsorbedSWFlux (1:sNx,1:sNy)
93     _RL IcePenetSWFluxFrac (1:sNx,1:sNy)
94    
95     C local copies of global variables
96     _RL tsurfLoc (1:sNx,1:sNy)
97     _RL atempLoc (1:sNx,1:sNy)
98     _RL lwdownLoc (1:sNx,1:sNy)
99     _RL ALB (1:sNx,1:sNy)
100     _RL ALB_ICE (1:sNx,1:sNy)
101     _RL ALB_SNOW (1:sNx,1:sNy)
102    
103 jmc 1.3 C i, j :: Loop counters
104     C kSrf :: vertical index of surface layer
105 jmc 1.1 INTEGER i, j
106 jmc 1.3 #ifdef SEAICE_VARIABLE_FREEZING_POINT
107     INTEGER kSrf
108     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
109 jmc 1.1 INTEGER ITER
110    
111 mlosch 1.5 C This is HICE_ACTUAL.GT.0.
112     LOGICAL iceOrNot(1:sNx,1:sNy)
113    
114     C TB :: temperature in boundary layer (=freezing point temperature) (K)
115     _RL TB (1:sNx,1:sNy)
116     C
117     _RL D1, D1I, D3
118     _RL TMELT, XKI, XKS, HCUT, XIO
119     _RL SurfMeltTemp
120     C effective conductivity of combined ice and snow
121     _RL effConduct(1:sNx,1:sNy)
122 jmc 1.1
123 jmc 1.3 C Constants to calculate Saturation Vapor Pressure
124 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
125 mlosch 1.5 _RL TMELTP, C1, C2, C3, C4, C5, QS1
126     _RL A2 (1:sNx,1:sNy)
127     _RL A3 (1:sNx,1:sNy)
128     c _RL B (1:sNx,1:sNy)
129     _RL A1 (1:sNx,1:sNy)
130 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
131 jmc 1.3 _RL dFiDTs1
132 jmc 1.1 _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
133     C specific humidity at ice surface variables
134 mlosch 1.5 _RL mm_pi,mm_log10pi,dqhice_dTice
135 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
136 jmc 1.1
137 ifenty 1.14 C latent heat of sublimation for ice (SEAICE_lhEvap +
138     C SEAICE_lhFusion)
139     _RL lhSublim
140    
141 jmc 1.1 C powers of temperature
142     _RL t1, t2, t3, t4
143 mlosch 1.5 _RL lnTEN
144     CEOP
145 jmc 1.1
146 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
147     CADJ INIT comlev1_solve4temp = COMMON, sNx*sNy*NMAX_TICE
148     #endif /* ALLOW_AUTODIFF_TAMC */
149    
150 mlosch 1.5 lnTEN = log(10.0 _d 0)
151 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
152 jmc 1.1 C MAYKUTS CONSTANTS FOR SAT. VAP. PRESSURE TEMP. POLYNOMIAL
153 mlosch 1.5 C1= 2.7798202 _d -06
154     C2= -2.6913393 _d -03
155     C3= 0.97920849 _d +00
156     C4= -158.63779 _d +00
157     C5= 9653.1925 _d +00
158 jmc 1.1
159     QS1=0.622 _d +00/1013.0 _d +00
160    
161 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
162 jmc 1.1 aa1 = 2663.5 _d 0
163     aa2 = 12.537 _d 0
164     bb1 = 0.622 _d 0
165 mlosch 1.5 bb2 = 1.0 _d 0 - bb1
166 jmc 1.1 Ppascals = 100000. _d 0
167 mlosch 1.5 C cc0 = TEN ** aa2
168     cc0 = exp(aa2*lnTEN)
169     cc1 = cc0*aa1*bb1*Ppascals*lnTEN
170 jmc 1.1 cc2 = cc0*bb2
171 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
172 jmc 1.1
173 jmc 1.4 #ifdef SEAICE_VARIABLE_FREEZING_POINT
174     kSrf = 1
175     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
176 jmc 1.1
177     C SENSIBLE HEAT CONSTANT
178 mlosch 1.7 D1=SEAICE_dalton*SEAICE_cpAir*SEAICE_rhoAir
179 jmc 1.1
180     C ICE LATENT HEAT CONSTANT
181 ifenty 1.14 lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
182     D1I=SEAICE_dalton*lhSublim*SEAICE_rhoAir
183 jmc 1.1
184     C STEFAN BOLTZMAN CONSTANT TIMES 0.97 EMISSIVITY
185     D3=SEAICE_emissivity
186    
187     C MELTING TEMPERATURE OF ICE
188 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
189 mlosch 1.5 TMELT = 273.16 _d +00
190     TMELTP = 273.159 _d +00
191 jmc 1.1 SurfMeltTemp = TMELTP
192 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
193 mlosch 1.5 TMELT = celsius2K
194 jmc 1.1 SurfMeltTemp = TMELT
195 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
196 jmc 1.1
197     C ICE CONDUCTIVITY
198     XKI=SEAICE_iceConduct
199    
200     C SNOW CONDUCTIVITY
201     XKS=SEAICE_snowConduct
202    
203     C CUTOFF SNOW THICKNESS
204     HCUT=SEAICE_snowThick
205    
206     C PENETRATION SHORTWAVE RADIATION FACTOR
207     XIO=SEAICE_shortwave
208    
209 jmc 1.3 C Initialize variables
210 jmc 1.1 DO J=1,sNy
211 mlosch 1.5 DO I=1,sNx
212     C HICE_ACTUAL is modified in this routine, but at the same time
213     C used to decided where there is ice, therefore we save this information
214     C here in a separate array
215     iceOrNot (I,J) = HICE_ACTUAL(I,J) .GT. 0. _d 0
216 jmc 1.6 C
217 mlosch 1.5 IcePenetSWFlux (I,J) = 0. _d 0
218     IcePenetSWFluxFrac (I,J) = 0. _d 0
219     AbsorbedSWFlux (I,J) = 0. _d 0
220    
221     qhice (I,J) = 0. _d 0
222     F_ia (I,J) = 0. _d 0
223 jmc 1.6
224 mlosch 1.5 F_io_net (I,J) = 0. _d 0
225     F_ia_net (I,J) = 0. _d 0
226    
227 mlosch 1.10 F_lh (I,J) = 0. _d 0
228    
229 mlosch 1.5 C Reset the snow/ice surface to TMELT and bound the atmospheric temperature
230 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
231 mlosch 1.5 tsurfLoc (I,J) = MIN(273.16 _d 0 + MAX_TICE,TSURF(I,J,bi,bj))
232     atempLoc (I,J) = MAX(273.16 _d 0 + MIN_ATEMP,ATEMP(I,J,bi,bj))
233     A1(I,J) = 0.0 _d 0
234     A2(I,J) = 0.0 _d 0
235     A3(I,J) = 0.0 _d 0
236     c B(I,J) = 0.0 _d 0
237     lwdownLoc(I,J) = MAX(MIN_LWDOWN,LWDOWN(I,J,bi,bj))
238 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
239 mlosch 1.5 F_swi (I,J) = 0. _d 0
240     F_lwd (I,J) = 0. _d 0
241     F_lwu (I,J) = 0. _d 0
242     F_sens (I,J) = 0. _d 0
243 jmc 1.6
244 mlosch 1.5 tsurfLoc (I,J) = TSURF(I,J,bi,bj)
245     atempLoc (I,J) = MAX(TMELT + MIN_ATEMP,ATEMP(I,J,bi,bj))
246     lwdownLoc(I,J) = LWDOWN(I,J,bi,bj)
247 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
248 jmc 1.1
249 mlosch 1.5 C FREEZING TEMPERATURE OF SEAWATER
250     #ifdef SEAICE_VARIABLE_FREEZING_POINT
251     C Use a variable seawater freezing point
252     TB(I,J) = -0.0575 _d 0*salt(I,J,kSrf,bi,bj) + 0.0901 _d 0
253     & + celsius2K
254     #else
255     C Use a constant freezing temperature (SEAICE_VARIABLE_FREEZING_POINT undef)
256 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
257 mlosch 1.5 TB(I,J) = 271.2 _d 0
258 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
259 mlosch 1.5 TB(I,J) = celsius2K + SEAICE_freeze
260 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
261 mlosch 1.5 #endif /* SEAICE_VARIABLE_FREEZING_POINT */
262     ENDDO
263 jmc 1.1 ENDDO
264    
265     DO J=1,sNy
266 mlosch 1.5 DO I=1,sNx
267 jmc 1.1
268     C DECIDE ON ALBEDO
269 mlosch 1.5 IF ( iceOrNot(I,J) ) THEN
270 jmc 1.6
271 mlosch 1.5 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 ) THEN
272     IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
273     ALB_ICE (I,J) = SEAICE_wetIceAlb_south
274     ALB_SNOW(I,J) = SEAICE_wetSnowAlb_south
275     ELSE ! no surface melting
276     ALB_ICE (I,J) = SEAICE_dryIceAlb_south
277     ALB_SNOW(I,J) = SEAICE_drySnowAlb_south
278     ENDIF
279     ELSE !/ Northern Hemisphere
280     IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
281     ALB_ICE (I,J) = SEAICE_wetIceAlb
282     ALB_SNOW(I,J) = SEAICE_wetSnowAlb
283     ELSE ! no surface melting
284     ALB_ICE (I,J) = SEAICE_dryIceAlb
285     ALB_SNOW(I,J) = SEAICE_drySnowAlb
286     ENDIF
287     ENDIF !/ Albedo for snow and ice
288    
289 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
290 mlosch 1.5 C If actual snow thickness exceeds the cutoff thickness, use the
291     C snow albedo
292     IF (HSNOW_ACTUAL(I,J) .GT. HCUT) THEN
293     ALB(I,J) = ALB_SNOW(I,J)
294    
295     C otherwise, use some combination of ice and snow albedo
296     C (What is the source of this formulation ?)
297     ELSE
298     ALB(I,J) = MIN(ALB_ICE(I,J) + HSNOW_ACTUAL(I,J)/HCUT*
299     & (ALB_SNOW(I,J) -ALB_ICE(I,J)),
300     & ALB_SNOW(I,J))
301     ENDIF
302    
303 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
304 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
305     ALB(I,J) = ALB_SNOW(I,J)
306     ELSE
307     ALB(I,J) = ALB_ICE(I,J)
308     ENDIF
309 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
310 jmc 1.1
311 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
312 mlosch 1.5 C NOW DETERMINE FIXED FORCING TERM IN HEAT BUDGET
313 jmc 1.1
314     #ifdef ALLOW_DOWNWARD_RADIATION
315 mlosch 1.5 IF(HSNOW_ACTUAL(I,J).GT.0.0) THEN
316     C NO SW PENETRATION WITH SNOW
317     A1(I,J)=(1.0 _d 0 - ALB(I,J))*SWDOWN(I,J,bi,bj)
318     & +lwdownLoc(I,J)*0.97 _d 0
319     & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
320     ELSE
321 jmc 1.1 C SW PENETRATION UNDER ICE
322 mlosch 1.5 A1(I,J)=(1.0 _d 0 - ALB(I,J))*SWDOWN(I,J,bi,bj)
323     & *(1.0 _d 0 - XIO*EXP(-1.5 _d 0*HICE_ACTUAL(I,J)))
324     & +lwdownLoc(I,J)*0.97 _d 0
325     & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
326     ENDIF
327 jmc 1.2 #endif /* ALLOW_DOWNWARD_RADIATION */
328 jmc 1.1
329 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
330 jmc 1.1
331 jmc 1.3 C The longwave radiative flux convergence
332 mlosch 1.5 F_lwd(I,J) = - 0.97 _d 0 * lwdownLoc(I,J)
333 jmc 1.1
334 jmc 1.3 C Determine the fraction of shortwave radiative flux
335     C remaining after scattering through the snow and ice at
336     C the ocean interface. If snow is present, no radiation
337     C penetrates to the ocean.
338 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
339     IcePenetSWFluxFrac(I,J) = 0.0 _d 0
340     ELSE
341     IcePenetSWFluxFrac(I,J) =
342     & XIO*EXP(-1.5 _d 0 * HICE_ACTUAL(I,J))
343     ENDIF
344 jmc 1.1
345 jmc 1.3 C The shortwave radiative flux convergence in the
346     C seaice.
347 mlosch 1.5 AbsorbedSWFlux(I,J) = -(1.0 _d 0 - ALB(I,J))*
348     & (1.0 _d 0 - IcePenetSWFluxFrac(I,J))
349     & *SWDOWN(I,J,bi,bj)
350 jmc 1.6
351 jmc 1.3 C The shortwave radiative flux convergence in the
352     C ocean beneath ice.
353 mlosch 1.5 IcePenetSWFlux(I,J) = -(1.0 _d 0 - ALB(I,J))*
354     & IcePenetSWFluxFrac(I,J)
355     & *SWDOWN(I,J,bi,bj)
356 jmc 1.1
357 mlosch 1.5 F_swi(I,J) = AbsorbedSWFlux(I,J)
358 jmc 1.6
359 jmc 1.3 C Set a mininum sea ice thickness of 5 cm to bound
360     C the magnitude of conductive heat fluxes.
361 ifenty 1.15 cif * now taken care of by SEAICE_hice_reg in seaice_growth
362     C hice_tmp = max(HICE_ACTUAL(I,J),5. _d -2)
363 jmc 1.1
364 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
365 jmc 1.1
366 jmc 1.3 C The effective conductivity of the two-layer
367     C snow/ice system.
368 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
369 mlosch 1.5 effConduct(I,J)=
370     & XKS/(HSNOW_ACTUAL(I,J)/HICE_ACTUAL(I,J) +
371     & XKS/XKI)/HICE_ACTUAL(I,J)
372 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
373 mlosch 1.5 effConduct(I,J) = XKI * XKS /
374 ifenty 1.15 & (XKS * HICE_ACTUAL(I,j) + XKI * HSNOW_ACTUAL(I,J))
375 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
376 jmc 1.1
377     #ifdef SEAICE_DEBUG
378 mlosch 1.5 IF ( (I .EQ. SEAICE_debugPointX) .and.
379     & (J .EQ. SEAICE_debugPointY) ) THEN
380 jmc 1.1
381 mlosch 1.5 print '(A,i6)','-----------------------------------'
382     print '(A,i6)','ibi merged initialization ', myIter
383 jmc 1.6
384 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
385     & 'ibi iter, TSL, TS ',myIter,
386     & tsurfLoc(I,J), TSURF(I,J,bi,bj)
387 jmc 1.6
388 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
389     & 'ibi iter, TMELT ',myIter,TMELT
390 jmc 1.6
391 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
392     & 'ibi iter, HIA, EFKCON ',myIter,
393     & HICE_ACTUAL(I,J), effConduct(I,J)
394 jmc 1.6
395 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
396     & 'ibi iter, HSNOW ',myIter,
397     & HSNOW_ACTUAL(I,J), ALB(I,J)
398 jmc 1.6
399 mlosch 1.5 print '(A,i6)','-----------------------------------'
400     print '(A,i6)','ibi energy balance iterat ', myIter
401 jmc 1.6
402 mlosch 1.5 ENDIF
403 jmc 1.2 #endif /* SEAICE_DEBUG */
404 jmc 1.6
405 mlosch 1.5 ENDIF !/* iceOrNot */
406     ENDDO !/* i */
407     ENDDO !/* j */
408 jmc 1.3 Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
409 mlosch 1.5 DO ITER=1,IMAX_TICE
410     DO J=1,sNy
411     DO I=1,sNx
412 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
413     iicekey = I + sNx*(J-1) + (ITER-1)*sNx*sNy
414     CADJ STORE tsurfloc(i,j) = comlev1_solve4temp,
415     CADJ & key = iicekey, byte = isbyte
416     #endif /* ALLOW_AUTODIFF_TAMC */
417    
418 mlosch 1.5 IF ( iceOrNot(I,J) ) THEN
419 jmc 1.1
420 mlosch 1.5 t1 = tsurfLoc(I,J)
421     t2 = t1*t1
422     t3 = t2*t1
423     t4 = t2*t2
424 jmc 1.1
425 mlosch 1.5 C Calculate the specific humidity in the BL above the snow/ice
426 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
427 mlosch 1.5 C Use the Maykut polynomial
428     qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
429 jmc 1.1
430 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
431 mlosch 1.5 C Use an approximation which is more accurate at low temperatures
432 jmc 1.1
433 mlosch 1.5 C log 10 of the sat vap pressure
434     mm_log10pi = -aa1 / t1 + aa2
435 jmc 1.1
436 mlosch 1.5 C The saturation vapor pressure (SVP) in the surface
437     C boundary layer (BL) above the snow/ice.
438     C mm_pi = TEN **(mm_log10pi)
439 jmc 1.6 C The following form does the same, but is faster
440 mlosch 1.5 mm_pi = exp(mm_log10pi*lnTEN)
441 jmc 1.1
442 mlosch 1.5 qhice(I,J) = bb1*mm_pi / (Ppascals - (1.0 _d 0 - bb1) *
443     & mm_pi)
444 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
445 jmc 1.1
446 mlosch 1.10 C Calculate the flux terms based on the updated tsurfLoc
447     F_c(I,J) = -effConduct(I,J)*(TB(I,J)-tsurfLoc(I,J))
448     F_lh(I,J) = D1I*UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
449 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
450 mlosch 1.5 A2(I,J)=-D1*UG(I,J)*t1-D1I*UG(I,J)*qhice(I,J)-D3*t4
451     A3(I,J) = 4.0 _d 0 * D3 * t3 + effConduct(I,J) + D1*UG(I,J)
452 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
453 mlosch 1.5 C A constant for SVP derivative w.r.t TICE
454     C cc3t = TEN **(aa1 / t1)
455 jmc 1.6 C The following form does the same, but is faster
456 mlosch 1.5 cc3t = exp(aa1 / t1 * lnTEN)
457    
458     c d(qh)/d(TICE)
459     dqhice_dTice = cc1*cc3t/((cc2-cc3t*Ppascals)**2 *t2)
460    
461     c d(F_ia)/d(TICE)
462     dFiDTs1 = 4.0 _d 0 * D3*t3 + effConduct(I,J) + D1*UG(I,J)
463     & + D1I*UG(I,J)*dqhice_dTice
464    
465     F_lwu(I,J)= t4 * D3
466 jmc 1.6
467 mlosch 1.5 F_sens(I,J)= D1 * UG(I,J) * (t1 - atempLoc(I,J))
468 jmc 1.6
469 mlosch 1.5 F_ia(I,J) = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +
470     & F_c(I,J) + F_sens(I,J) + F_lh(I,J)
471 jmc 1.1
472 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
473 jmc 1.1
474     #ifdef SEAICE_DEBUG
475 mlosch 1.5 IF ( (I .EQ. SEAICE_debugPointX) .and.
476     & (J .EQ. SEAICE_debugPointY) ) THEN
477     print '(A,i6,4(1x,D24.15))',
478     & 'ice-iter qhICE, ', ITER,qhIce(I,J)
479 jmc 1.1
480 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
481 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
482     & 'ice-iter A1 A2 B ', ITER,A1(I,J), A2(I,J),
483     & -F_c(I,J)
484 jmc 1.1
485 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
486     & 'ice-iter A3 (-A1+A2) ', ITER, A3(I,J),
487     & -(A1(I,J) + A2(I,J))
488 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
489    
490 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
491     & 'ice-iter dFiDTs1 F_ia ', ITER, dFiDTs1,
492     & F_ia(I,J)
493 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
494 jmc 1.1
495 mlosch 1.5 ENDIF
496 jmc 1.2 #endif /* SEAICE_DEBUG */
497 jmc 1.1
498 mlosch 1.5 C Update tsurfLoc
499 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
500 mlosch 1.5 tsurfLoc(I,J)=tsurfLoc(I,J)
501     & +(A1(I,J)+A2(I,J)-F_c(I,J))/A3(I,J)
502 jmc 1.6
503 mlosch 1.5 tsurfLoc(I,J) =MAX(273.16 _d 0+MIN_TICE,tsurfLoc(I,J))
504     tsurfLoc(I,J) =MIN(tsurfLoc(I,J),TMELT)
505 jmc 1.6
506     #else /* SEAICE_SOLVE4TEMP_LEGACY */
507 mlosch 1.5 tsurfLoc(I,J) = tsurfLoc(I,J) - F_ia(I,J) / dFiDTs1
508 jmc 1.1
509 jmc 1.3 C If the search leads to tsurfLoc < 50 Kelvin,
510     C restart the search at tsurfLoc = TMELT. Note that one
511     C solution to the energy balance problem is an
512     C extremely low temperature - a temperature far below
513     C realistic values.
514 jmc 1.1
515 mlosch 1.5 IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) THEN
516     tsurfLoc(I,J) = TMELT
517     ENDIF
518 gforget 1.11 tsurfLoc(I,J) =MIN(tsurfLoc(I,J),TMELT)
519 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
520 jmc 1.1
521     #ifdef SEAICE_DEBUG
522 mlosch 1.5 IF ( (I .EQ. SEAICE_debugPointX) .and.
523     & (J .EQ. SEAICE_debugPointY) ) THEN
524 jmc 1.6
525 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
526     & 'ice-iter tsurfLc,|dif|', ITER,
527     & tsurfLoc(I,J),
528     & log10(abs(tsurfLoc(I,J) - t1))
529     ENDIF
530 jmc 1.2 #endif /* SEAICE_DEBUG */
531 jmc 1.1
532 mlosch 1.5 ENDIF !/* iceOrNot */
533     ENDDO !/* i */
534     ENDDO !/* j */
535     ENDDO !/* Iterations */
536 jmc 1.3 Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
537 mlosch 1.5 DO J=1,sNy
538     DO I=1,sNx
539     IF ( iceOrNot(I,J) ) THEN
540 jmc 1.1
541 jmc 1.3 C Finalize the flux terms
542 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
543 mlosch 1.5 F_ia(I,J)=-A1(I,J)-A2(I,J)
544     TSURF(I,J,bi,bj)=MIN(tsurfLoc(I,J),TMELT)
545 jmc 1.6
546 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0 ) THEN
547     C NO SW PENETRATION WITH SNOW
548     IcePenetSWFlux(I,J)=0.0 _d 0
549     ELSE
550     C SW PENETRATION UNDER ICE
551 jmc 1.1
552     #ifdef ALLOW_DOWNWARD_RADIATION
553 mlosch 1.5 IcePenetSWFlux(I,J)=-(1.0 _d 0 -ALB(I,J))*SWDOWN(I,J,bi,bj)
554     & *XIO*EXP(-1.5 _d 0*HICE_ACTUAL(I,J))
555 jmc 1.2 #endif /* ALLOW_DOWNWARD_RADIATION */
556 mlosch 1.5 ENDIF
557 jmc 1.1
558 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
559 mlosch 1.5 TSURF(I,J,bi,bj) = tsurfLoc(I,J)
560 jmc 1.1
561 mlosch 1.5 C Recalculate the fluxes based on the (possibly) adjusted TSURF
562     t1 = tsurfLoc(I,J)
563     t2 = t1*t1
564     t3 = t2*t1
565     t4 = t2*t2
566 jmc 1.1
567 mlosch 1.5 C log 10 of the sat vap pressure
568     mm_log10pi = -aa1 / t1 + aa2
569 jmc 1.1
570 mlosch 1.5 C saturation vapor pressure
571     C mm_pi = TEN **(mm_log10pi)
572 jmc 1.6 C The following form does the same, but is faster
573 mlosch 1.5 mm_pi = exp(mm_log10pi*lnTEN)
574 jmc 1.1
575 jmc 1.3 C over ice specific humidity
576 mlosch 1.5 qhice(I,J) = bb1*mm_pi/(Ppascals- (1.0 _d 0 - bb1) * mm_pi)
577 jmc 1.1
578 mlosch 1.5 F_lh(I,J) = D1I * UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
579     F_c(I,J) = -effConduct(I,J) * (TB(I,J) - t1)
580     F_lwu(I,J) = t4 * D3
581     F_sens(I,J) = D1 * UG(I,J) * (t1 - atempLoc(I,J))
582 jmc 1.1
583 jmc 1.3 C The flux between the ice/snow surface and the atmosphere.
584     C (excludes upward conductive fluxes)
585 mlosch 1.5 F_ia(I,J) = F_lwd(I,J) + F_swi(I,J) + F_lwu(I,J) +
586     & F_sens(I,J) + F_lh(I,J)
587 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
588 jmc 1.1
589 gforget 1.9 #ifdef SEAICE_MODIFY_GROWTH_ADJ
590     Cgf no additional dependency through solver, snow, etc.
591     if ( SEAICEadjMODE.GE.2 ) then
592     CALL ZERO_ADJ_1D( 1, TSURF(I,J,bi,bj), myThid)
593     t1 = TSURF(I,J,bi,bj)
594     t2 = t1*t1
595     t3 = t2*t1
596     t4 = t2*t2
597     qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
598    
599     A1(I,J)=0.3 _d 0 *SWDOWN(I,J,bi,bj)+lwdownLoc(I,J)*0.97 _d 0
600     & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
601     A2(I,J)=-D1*UG(I,J)*t1-D1I*UG(I,J)*qhice(I,J)-D3*t4
602    
603     F_ia(I,J)=-A1(I,J)-A2(I,J)
604     IcePenetSWFlux(I,J)= 0. _d 0
605     endif
606     #endif
607    
608 mlosch 1.5 C Caclulate the net ice-ocean and ice-atmosphere fluxes
609     IF (F_c(I,J) .LT. 0.0 _d 0) THEN
610     F_io_net(I,J) = -F_c(I,J)
611     F_ia_net(I,J) = 0.0 _d 0
612     ELSE
613     F_io_net(I,J) = 0.0 _d 0
614     F_ia_net(I,J) = F_ia(I,J)
615     ENDIF !/* conductive fluxes up or down */
616 mlosch 1.10 C Fresh water flux (kg/m^2/s) from latent heat of sublimation.
617     C F_lh is positive upward (sea ice looses heat) and FWsublim
618     C is also positive upward (atmosphere gains freshwater)
619 ifenty 1.14 FWsublim(I,J) = F_lh(I,J)/lhSublim
620 jmc 1.1
621     #ifdef SEAICE_DEBUG
622 mlosch 1.5 IF ( (I .EQ. SEAICE_debugPointX) .and.
623     & (J .EQ. SEAICE_debugPointY) ) THEN
624 jmc 1.6
625 mlosch 1.5 print '(A)','----------------------------------------'
626     print '(A,i6)','ibi complete ', myIter
627 jmc 1.6
628 mlosch 1.5 print '(A,4(1x,D24.15))',
629     & 'ibi T(SURF, surfLoc,atmos) ',
630     & TSURF(I,J,bi,bj), tsurfLoc(I,J),atempLoc(I,J)
631 jmc 1.6
632 mlosch 1.5 print '(A,4(1x,D24.15))',
633     & 'ibi LWL ', lwdownLoc(I,J)
634 jmc 1.6
635 mlosch 1.5 print '(A,4(1x,D24.15))',
636     & 'ibi QSW(Total, Penetrating)',
637     & SWDOWN(I,J,bi,bj), IcePenetSWFlux(I,J)
638 jmc 1.6
639 mlosch 1.5 print '(A,4(1x,D24.15))',
640     & 'ibi qh(ATM ICE) ',
641     & AQH(I,J,bi,bj),qhice(I,J)
642 jmc 1.6
643 mlosch 1.5 c print '(A,4(1x,D24.15))',
644     c & 'ibi F(lwd,swi,lwu) ',
645     c & F_lwd(I,J), F_swi(I,J), F_lwu(I,J)
646    
647     c print '(A,4(1x,D24.15))',
648     c & 'ibi F(c,lh,sens) ',
649     c & F_c(I,J), F_lh(I,J), F_sens(I,J)
650    
651     print '(A,4(1x,D24.15))',
652     & 'ibi F_ia, F_ia_net, F_c ',
653 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
654 mlosch 1.5 & -(A1(I,J)+A2(I,J)),
655     & -(A1(I,J)+A2(I,J)-F_c(I,J)),
656     & F_c(I,J)
657 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
658 mlosch 1.5 & F_ia(I,J),
659     & F_ia_net(I,J),
660     & F_c(I,J)
661 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
662 jmc 1.1
663 mlosch 1.5 print '(A)','----------------------------------------'
664 jmc 1.1
665 mlosch 1.5 ENDIF
666 jmc 1.2 #endif /* SEAICE_DEBUG */
667 jmc 1.6
668 mlosch 1.5 ENDIF !/* iceOrNot */
669     ENDDO !/* i */
670 jmc 1.1 ENDDO !/* j */
671    
672     RETURN
673     END
  ViewVC Help
Powered by ViewVC 1.1.22