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