/[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.20 - (hide annotations) (download)
Wed Dec 28 20:54:28 2011 UTC (12 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63h, checkpoint63i
Changes since 1.19: +1 -2 lines
remove unused variables

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

  ViewVC Help
Powered by ViewVC 1.1.22