/[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.23 - (hide annotations) (download)
Tue Jan 31 15:57:17 2012 UTC (12 years, 4 months ago) by jmc
Branch: MAIN
Changes since 1.22: +11 -12 lines 
ifdef SOLVE4TEMP_LEGACY part:
- replace hard coded 273.16 (=previous value of celsius2K) by celsius2K
- replace hard coded 271.20 (=previous value of celsius2K + default
  value of SEAICE_freeze) by celsius2K+SEAICE_freeze so that seawater
  freezing temp is consistent with value in seaice_growth.F
1 jmc 1.23 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_solve4temp.F,v 1.22 2012/01/30 20:51:43 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 jmc 1.21 O F_ia, IcePenetSW,
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 jmc 1.21 C !INPUT PARAMETERS:
53     C UG :: atmospheric wind speed (m/s)
54 jmc 1.1 C HICE_ACTUAL :: actual ice thickness
55     C HSNOW_ACTUAL :: actual snow thickness
56 jmc 1.21 C TSURF :: surface temperature of ice/snow in Kelvin
57     C bi,bj :: tile indices
58     C myTime :: current time in simulation
59     C myIter :: iteration number in simulation
60     C myThid :: my Thread Id number
61     C !OUTPUT PARAMETERS:
62     C TSURF :: updated surface temperature of ice/snow in Kelvin
63     C F_ia :: upward seaice/snow surface heat flux to atmosphere (W/m^2)
64     C IcePenetSW :: short wave heat flux transmitted through ice (+=upward)
65     C FWsublim :: fresh water (mass) flux due to sublimation (+=up)(kg/m^2/s)
66     _RL UG (1:sNx,1:sNy)
67     _RL HICE_ACTUAL (1:sNx,1:sNy)
68     _RL HSNOW_ACTUAL(1:sNx,1:sNy)
69 ifenty 1.16 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
70 jmc 1.21 _RL F_lh_max (1:sNx,1:sNy)
71 ifenty 1.16 #endif
72 jmc 1.21 _RL TSURF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
73     _RL F_ia (1:sNx,1:sNy)
74     _RL IcePenetSW (1:sNx,1:sNy)
75     _RL FWsublim (1:sNx,1:sNy)
76 jmc 1.1 INTEGER bi, bj
77     _RL myTime
78     INTEGER myIter, myThid
79 jmc 1.19 CEOP
80 jmc 1.1
81 jmc 1.19 #if defined(ALLOW_ATM_TEMP) && defined(ALLOW_DOWNWARD_RADIATION)
82 jmc 1.1 C !LOCAL VARIABLES:
83     C === Local variables ===
84 jmc 1.3 C i, j :: Loop counters
85     C kSrf :: vertical index of surface layer
86 jmc 1.1 INTEGER i, j
87 jmc 1.3 #ifdef SEAICE_VARIABLE_FREEZING_POINT
88     INTEGER kSrf
89     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
90 jmc 1.1 INTEGER ITER
91 jmc 1.21 C TB :: ocean temperature in contact with ice (=seawater freezing point) (K)
92 mlosch 1.5 _RL TB (1:sNx,1:sNy)
93 mlosch 1.18 _RL D1, D1I
94     _RL D3(1:sNx,1:sNy)
95 mlosch 1.5 _RL TMELT, XKI, XKS, HCUT, XIO
96     _RL SurfMeltTemp
97 jmc 1.21 C effConduct :: effective conductivity of combined ice and snow
98 mlosch 1.5 _RL effConduct(1:sNx,1:sNy)
99 jmc 1.21 C lhSublim :: latent heat of sublimation (SEAICE_lhEvap + SEAICE_lhFusion)
100     _RL lhSublim
101     C t1,t2,t3,t4 :: powers of temperature
102     _RL t1, t2, t3, t4
103 jmc 1.1
104 jmc 1.3 C Constants to calculate Saturation Vapor Pressure
105 jmc 1.21 _RL dqhice_dTice
106 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
107 jmc 1.21 _RL C1, C2, C3, C4, C5, QS1
108 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
109 jmc 1.3 _RL dFiDTs1
110 jmc 1.21 _RL lnTEN
111 jmc 1.1 _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
112     C specific humidity at ice surface variables
113 jmc 1.21 _RL mm_pi,mm_log10pi
114 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
115 jmc 1.1
116 jmc 1.22 C F_c :: conductive heat flux through seaice+snow (+=upward)
117 jmc 1.21 C F_lh :: latent heat flux (sublimation) (+=upward)
118     #ifdef SEAICE_SOLVE4TEMP_LEGACY
119     C A1 :: part of atmos surface flux (+=downward) independent of tsurf
120     C A2 :: part of atmos surface flux (+=upward) which depends on tsurf
121 jmc 1.22 C A3 :: derivative of (A2-F_c) versus tsurf
122 jmc 1.21 _RL A2 (1:sNx,1:sNy)
123     _RL A3 (1:sNx,1:sNy)
124     _RL A1 (1:sNx,1:sNy)
125     #else /* SEAICE_SOLVE4TEMP_LEGACY */
126 jmc 1.22 C F_lwu :: upward long-wave surface heat flux (+=upward)
127     C F_sens :: sensible surface heat flux (+=upward)
128 jmc 1.21 _RL F_lwu (1:sNx,1:sNy)
129     _RL F_sens (1:sNx,1:sNy)
130     #endif /* SEAICE_SOLVE4TEMP_LEGACY */
131     _RL F_lh (1:sNx,1:sNy)
132     _RL F_c (1:sNx,1:sNy)
133     _RL qhice (1:sNx,1:sNy)
134     _RL absorbedSW (1:sNx,1:sNy)
135     _RL penetSWFrac
136    
137     C local copies of global variables
138     _RL tsurfLoc (1:sNx,1:sNy)
139     _RL atempLoc (1:sNx,1:sNy)
140     _RL lwdownLoc (1:sNx,1:sNy)
141     _RL ALB (1:sNx,1:sNy)
142     _RL ALB_ICE (1:sNx,1:sNy)
143     _RL ALB_SNOW (1:sNx,1:sNy)
144     C iceOrNot :: this is HICE_ACTUAL.GT.0.
145     LOGICAL iceOrNot(1:sNx,1:sNy)
146     #ifdef SEAICE_DEBUG
147     C F_io_net :: upward conductive heat flux through seaice+snow
148     C F_ia_net :: net heat flux divergence at the sea ice/snow surface:
149     C includes ice conductive fluxes and atmospheric fluxes (W/m^2)
150     _RL F_io_net (1:sNx,1:sNy)
151     _RL F_ia_net (1:sNx,1:sNy)
152     #endif /* SEAICE_DEBUG */
153 ifenty 1.14
154 jmc 1.21 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
155 jmc 1.1
156 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
157     CADJ INIT comlev1_solve4temp = COMMON, sNx*sNy*NMAX_TICE
158     #endif /* ALLOW_AUTODIFF_TAMC */
159    
160 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
161 jmc 1.21 C MAYKUT 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 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
170 jmc 1.21 lnTEN = LOG(10.0 _d 0)
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 jmc 1.21 cc0 = EXP(aa2*lnTEN)
178 mlosch 1.5 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 mlosch 1.5 TMELT = celsius2K
195 jmc 1.1 SurfMeltTemp = TMELT
196 jmc 1.23 #ifdef SEAICE_SOLVE4TEMP_LEGACY
197     C old SOLVE4TEMP_LEGACY setting, consistent with former celsius2K value:
198     c TMELT = 273.16 _d +00
199     c SurfMeltTemp = 273.159 _d +00
200     SurfMeltTemp = TMELT - 1. _d -3
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 jmc 1.21 C Snow-Thickness above HCUT: SW optically thick snow (=> snow-albedo).
211     C Snow-Thickness below HCUT: linear transition to ice-albedo
212     #ifdef SEAICE_SOLVE4TEMP_LEGACY
213     HCUT = SEAICE_snowThick
214     #else /* SEAICE_SOLVE4TEMP_LEGACY */
215     HCUT = 0. _d 0
216     #endif /* SEAICE_SOLVE4TEMP_LEGACY */
217 jmc 1.1
218     C PENETRATION SHORTWAVE RADIATION FACTOR
219     XIO=SEAICE_shortwave
220    
221 jmc 1.3 C Initialize variables
222 jmc 1.1 DO J=1,sNy
223 mlosch 1.5 DO I=1,sNx
224     C HICE_ACTUAL is modified in this routine, but at the same time
225     C used to decided where there is ice, therefore we save this information
226     C here in a separate array
227 jmc 1.21 iceOrNot (I,J) = HICE_ACTUAL(I,J) .GT. 0. _d 0
228     IcePenetSW(I,J) = 0. _d 0
229     absorbedSW(I,J) = 0. _d 0
230 mlosch 1.5 qhice (I,J) = 0. _d 0
231     F_ia (I,J) = 0. _d 0
232 jmc 1.21 c F_io_net (I,J) = 0. _d 0
233     c F_ia_net (I,J) = 0. _d 0
234 mlosch 1.10 F_lh (I,J) = 0. _d 0
235    
236 mlosch 1.5 C Reset the snow/ice surface to TMELT and bound the atmospheric temperature
237 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
238 mlosch 1.5 A1(I,J) = 0.0 _d 0
239     A2(I,J) = 0.0 _d 0
240     A3(I,J) = 0.0 _d 0
241 jmc 1.23 tsurfLoc (I,J) = MIN( celsius2K+MAX_TICE, TSURF(I,J,bi,bj) )
242     atempLoc (I,J) = MAX( celsius2K+MIN_ATEMP, ATEMP(I,J,bi,bj) )
243 jmc 1.21 lwdownLoc(I,J) = MAX( MIN_LWDOWN, LWDOWN(I,J,bi,bj) )
244 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
245 mlosch 1.5 F_lwu (I,J) = 0. _d 0
246     F_sens (I,J) = 0. _d 0
247     tsurfLoc (I,J) = TSURF(I,J,bi,bj)
248 jmc 1.21 atempLoc (I,J) = MAX( celsius2K+MIN_ATEMP, ATEMP(I,J,bi,bj) )
249 mlosch 1.5 lwdownLoc(I,J) = LWDOWN(I,J,bi,bj)
250 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
251 jmc 1.1
252 mlosch 1.5 C FREEZING TEMPERATURE OF SEAWATER
253     #ifdef SEAICE_VARIABLE_FREEZING_POINT
254     C Use a variable seawater freezing point
255     TB(I,J) = -0.0575 _d 0*salt(I,J,kSrf,bi,bj) + 0.0901 _d 0
256     & + celsius2K
257     #else
258     C Use a constant freezing temperature (SEAICE_VARIABLE_FREEZING_POINT undef)
259 jmc 1.23 C old SOLVE4TEMP_LEGACY setting (not consistent with seaice_growth value)
260     c TB(I,J) = 271.2 _d 0
261 mlosch 1.5 TB(I,J) = celsius2K + SEAICE_freeze
262     #endif /* SEAICE_VARIABLE_FREEZING_POINT */
263 mlosch 1.18 IF(HSNOW_ACTUAL(I,J).GT.0.0) THEN
264 jmc 1.21 C Stefan-Boltzmann constant times emissivity
265 mlosch 1.18 D3(I,J)=SEAICE_snow_emiss*SEAICE_boltzmann
266     #ifdef EXF_LWDOWN_WITH_EMISSIVITY
267     C This is now [(1-emiss)*lwdown - lwdown]
268     lwdownloc(I,J) = SEAICE_snow_emiss*lwdownloc(I,J)
269     #else /* use the old hard wired inconsistent value */
270     lwdownloc(I,J) = 0.97 _d 0*lwdownloc(I,J)
271     #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
272     ELSE
273 jmc 1.21 C Stefan-Boltzmann constant times emissivity
274 mlosch 1.18 D3(I,J)=SEAICE_ice_emiss*SEAICE_boltzmann
275     #ifdef EXF_LWDOWN_WITH_EMISSIVITY
276     C This is now [(1-emiss)*lwdown - lwdown]
277     lwdownloc(I,J) = SEAICE_ice_emiss*lwdownloc(I,J)
278     #else /* use the old hard wired inconsistent value */
279     lwdownloc(I,J) = 0.97 _d 0*lwdownloc(I,J)
280     #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
281     ENDIF
282 mlosch 1.5 ENDDO
283 jmc 1.1 ENDDO
284    
285     DO J=1,sNy
286 mlosch 1.5 DO I=1,sNx
287 jmc 1.1
288     C DECIDE ON ALBEDO
289 mlosch 1.5 IF ( iceOrNot(I,J) ) THEN
290 jmc 1.6
291 mlosch 1.5 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 ) THEN
292     IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
293     ALB_ICE (I,J) = SEAICE_wetIceAlb_south
294     ALB_SNOW(I,J) = SEAICE_wetSnowAlb_south
295     ELSE ! no surface melting
296     ALB_ICE (I,J) = SEAICE_dryIceAlb_south
297     ALB_SNOW(I,J) = SEAICE_drySnowAlb_south
298     ENDIF
299     ELSE !/ Northern Hemisphere
300     IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
301     ALB_ICE (I,J) = SEAICE_wetIceAlb
302     ALB_SNOW(I,J) = SEAICE_wetSnowAlb
303     ELSE ! no surface melting
304     ALB_ICE (I,J) = SEAICE_dryIceAlb
305     ALB_SNOW(I,J) = SEAICE_drySnowAlb
306     ENDIF
307     ENDIF !/ Albedo for snow and ice
308    
309 jmc 1.21 C If actual snow thickness exceeds the cutoff thickness, use snow albedo
310 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. HCUT) THEN
311     ALB(I,J) = ALB_SNOW(I,J)
312 jmc 1.21 ELSEIF ( HCUT.LE.ZERO ) THEN
313     ALB(I,J) = ALB_ICE(I,J)
314 mlosch 1.5 ELSE
315 jmc 1.21 C otherwise, use linear transition between ice and snow albedo
316     ALB(I,J) = MIN( ALB_ICE(I,J) + HSNOW_ACTUAL(I,J)/HCUT
317     & *(ALB_SNOW(I,J) -ALB_ICE(I,J))
318     & , ALB_SNOW(I,J) )
319 mlosch 1.5 ENDIF
320    
321 jmc 1.21 C Determine the fraction of shortwave radiative flux remaining
322     C at ocean interface after scattering through the snow and ice.
323     C If snow is present, no radiation penetrates through snow+ice
324 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
325 jmc 1.21 penetSWFrac = 0.0 _d 0
326 mlosch 1.5 ELSE
327 jmc 1.21 penetSWFrac = XIO*EXP(-1.5 _d 0 * HICE_ACTUAL(I,J))
328 mlosch 1.5 ENDIF
329 jmc 1.21 C The shortwave radiative flux leaving ocean beneath ice (+=up).
330     IcePenetSW(I,J) = -(1.0 _d 0 - ALB(I,J))
331     & *penetSWFrac * SWDOWN(I,J,bi,bj)
332     C The shortwave radiative flux convergence in the seaice.
333     absorbedSW(I,J) = (1.0 _d 0 - ALB(I,J))
334     & *(1.0 _d 0 - penetSWFrac)* SWDOWN(I,J,bi,bj)
335 jmc 1.1
336 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
337 jmc 1.21 C Now determine fixed (relative to tsurf) forcing term in heat budget
338     A1(I,J) = absorbedSW(I,J) + lwdownLoc(I,J)
339 mlosch 1.5 & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
340 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
341 jmc 1.1
342 jmc 1.3 C The effective conductivity of the two-layer
343     C snow/ice system.
344 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
345 mlosch 1.5 effConduct(I,J)=
346     & XKS/(HSNOW_ACTUAL(I,J)/HICE_ACTUAL(I,J) +
347     & XKS/XKI)/HICE_ACTUAL(I,J)
348 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
349 jmc 1.21 C Set a mininum sea ice thickness of 5 cm to bound
350     C the magnitude of conductive heat fluxes.
351     Cif * now taken care of by SEAICE_hice_reg in seaice_growth
352     c hice_tmp = max(HICE_ACTUAL(I,J),5. _d -2)
353 mlosch 1.5 effConduct(I,J) = XKI * XKS /
354 jmc 1.21 & (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
355 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
356 jmc 1.1
357     #ifdef SEAICE_DEBUG
358 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
359 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
360 mlosch 1.5 print '(A,i6)','-----------------------------------'
361     print '(A,i6)','ibi merged initialization ', myIter
362     print '(A,i6,4(1x,D24.15))',
363     & 'ibi iter, TSL, TS ',myIter,
364     & tsurfLoc(I,J), TSURF(I,J,bi,bj)
365     print '(A,i6,4(1x,D24.15))',
366     & 'ibi iter, TMELT ',myIter,TMELT
367     print '(A,i6,4(1x,D24.15))',
368     & 'ibi iter, HIA, EFKCON ',myIter,
369     & HICE_ACTUAL(I,J), effConduct(I,J)
370     print '(A,i6,4(1x,D24.15))',
371     & 'ibi iter, HSNOW ',myIter,
372     & HSNOW_ACTUAL(I,J), ALB(I,J)
373     print '(A,i6)','-----------------------------------'
374     print '(A,i6)','ibi energy balance iterat ', myIter
375     ENDIF
376 jmc 1.2 #endif /* SEAICE_DEBUG */
377 jmc 1.6
378 mlosch 1.5 ENDIF !/* iceOrNot */
379     ENDDO !/* i */
380     ENDDO !/* j */
381 jmc 1.21
382     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
383 mlosch 1.5 DO ITER=1,IMAX_TICE
384     DO J=1,sNy
385     DO I=1,sNx
386 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
387     iicekey = I + sNx*(J-1) + (ITER-1)*sNx*sNy
388     CADJ STORE tsurfloc(i,j) = comlev1_solve4temp,
389     CADJ & key = iicekey, byte = isbyte
390     #endif /* ALLOW_AUTODIFF_TAMC */
391    
392 mlosch 1.5 IF ( iceOrNot(I,J) ) THEN
393 jmc 1.1
394 mlosch 1.5 t1 = tsurfLoc(I,J)
395     t2 = t1*t1
396     t3 = t2*t1
397     t4 = t2*t2
398 jmc 1.1
399 mlosch 1.5 C Calculate the specific humidity in the BL above the snow/ice
400 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
401 mlosch 1.5 C Use the Maykut polynomial
402     qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
403 jmc 1.21 dqhice_dTice = 0. _d 0
404 jmc 1.1
405 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
406 mlosch 1.5 C Use an approximation which is more accurate at low temperatures
407 jmc 1.1
408 mlosch 1.5 C log 10 of the sat vap pressure
409     mm_log10pi = -aa1 / t1 + aa2
410     C The saturation vapor pressure (SVP) in the surface
411     C boundary layer (BL) above the snow/ice.
412     C mm_pi = TEN **(mm_log10pi)
413 jmc 1.6 C The following form does the same, but is faster
414 jmc 1.21 mm_pi = EXP(mm_log10pi*lnTEN)
415 mlosch 1.5 qhice(I,J) = bb1*mm_pi / (Ppascals - (1.0 _d 0 - bb1) *
416     & mm_pi)
417 jmc 1.21 C A constant for SVP derivative w.r.t TICE
418     C cc3t = TEN **(aa1 / t1)
419     C The following form does the same, but is faster
420     cc3t = EXP(aa1 / t1 * lnTEN)
421     C d(qh)/d(TICE)
422     dqhice_dTice = cc1*cc3t/((cc2-cc3t*Ppascals)**2 *t2)
423 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
424 jmc 1.1
425 mlosch 1.10 C Calculate the flux terms based on the updated tsurfLoc
426 jmc 1.22 F_c(I,J) = effConduct(I,J)*(TB(I,J)-tsurfLoc(I,J))
427 mlosch 1.10 F_lh(I,J) = D1I*UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
428 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
429 jmc 1.21 A2(I,J) = D1*UG(I,J)*t1+D1I*UG(I,J)*qhice(I,J)+D3(I,J)*t4
430 mlosch 1.18 A3(I,J) = 4.0 _d 0*D3(I,J)*t3 + effConduct(I,J)+D1*UG(I,J)
431 jmc 1.21 & + D1I*UG(I,J)*dqhice_dTice
432 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
433 ifenty 1.16 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
434 jmc 1.21 C if the latent heat flux implied by tsurfLoc exceeds
435     C F_lh_max, cap F_lh and decouple the flux magnitude from TICE
436     IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
437 ifenty 1.16 F_lh(I,J) = F_lh_max(I,J)
438     dqhice_dTice = ZERO
439 jmc 1.21 ENDIF
440     #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
441 ifenty 1.16
442 jmc 1.21 C d(F_ia)/d(TICE)
443 mlosch 1.18 dFiDTs1 = 4.0 _d 0*D3(I,J)*t3 + effConduct(I,J) + D1*UG(I,J)
444 mlosch 1.5 & + D1I*UG(I,J)*dqhice_dTice
445    
446 jmc 1.21 F_lwu(I,J) = t4 * D3(I,J)
447 mlosch 1.5 F_sens(I,J)= D1 * UG(I,J) * (t1 - atempLoc(I,J))
448 jmc 1.21 F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
449 jmc 1.22 & -F_c(I,J) + F_sens(I,J) + F_lh(I,J)
450 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
451 jmc 1.1
452     #ifdef SEAICE_DEBUG
453 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
454 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
455 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
456     & 'ice-iter qhICE, ', ITER,qhIce(I,J)
457 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
458 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
459 jmc 1.22 & 'ice-iter A1 A2 B ',ITER,A1(I,J),A2(I,J),F_c(I,J)
460 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
461 jmc 1.21 & 'ice-iter A3 (-A1+A2) ',ITER,A3(I,J),-A1(I,J)+A2(I,J)
462 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
463 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
464 jmc 1.21 & 'ice-iter dFiDTs1 F_ia ', ITER, dFiDTs1, F_ia(I,J)
465 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
466 mlosch 1.5 ENDIF
467 jmc 1.2 #endif /* SEAICE_DEBUG */
468 jmc 1.1
469 mlosch 1.5 C Update tsurfLoc
470 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
471 jmc 1.22 C update tsurf as solution of : Fc = A2 - A1 + A3 *delta.tsurf
472 mlosch 1.5 tsurfLoc(I,J)=tsurfLoc(I,J)
473 jmc 1.22 & +(A1(I,J)-A2(I,J)+F_c(I,J))/A3(I,J)
474 jmc 1.23 tsurfLoc(I,J) = MAX( celsius2K+MIN_TICE, tsurfLoc(I,J) )
475 jmc 1.6
476     #else /* SEAICE_SOLVE4TEMP_LEGACY */
477 mlosch 1.5 tsurfLoc(I,J) = tsurfLoc(I,J) - F_ia(I,J) / dFiDTs1
478 jmc 1.1
479 jmc 1.21 C If the search leads to tsurfLoc < 50 Kelvin, restart the search
480     C at tsurfLoc = TMELT. Note that one solution to the energy balance problem
481     C is an extremely low temperature - a temperature far below realistic values.
482 mlosch 1.5 IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) THEN
483     tsurfLoc(I,J) = TMELT
484     ENDIF
485 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
486 jmc 1.21 tsurfLoc(I,J) = MIN( tsurfLoc(I,J), TMELT )
487 jmc 1.1
488     #ifdef SEAICE_DEBUG
489 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
490 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
491 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
492     & 'ice-iter tsurfLc,|dif|', ITER,
493     & tsurfLoc(I,J),
494 jmc 1.21 & LOG10(ABS(tsurfLoc(I,J) - t1))
495 mlosch 1.5 ENDIF
496 jmc 1.2 #endif /* SEAICE_DEBUG */
497 jmc 1.1
498 mlosch 1.5 ENDIF !/* iceOrNot */
499     ENDDO !/* i */
500     ENDDO !/* j */
501     ENDDO !/* Iterations */
502 jmc 1.21 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
503    
504 mlosch 1.5 DO J=1,sNy
505     DO I=1,sNx
506     IF ( iceOrNot(I,J) ) THEN
507 jmc 1.1
508 jmc 1.21 C Save updated tsurf and finalize the flux terms
509     TSURF(I,J,bi,bj) = tsurfLoc(I,J)
510    
511 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
512 jmc 1.21 F_ia(I,J)=-A1(I,J)+A2(I,J)
513 jmc 1.6 #else /* SEAICE_SOLVE4TEMP_LEGACY */
514 mlosch 1.5 C Recalculate the fluxes based on the (possibly) adjusted TSURF
515     t1 = tsurfLoc(I,J)
516     t2 = t1*t1
517     t3 = t2*t1
518     t4 = t2*t2
519 jmc 1.1
520 mlosch 1.5 C log 10 of the sat vap pressure
521     mm_log10pi = -aa1 / t1 + aa2
522     C saturation vapor pressure
523     C mm_pi = TEN **(mm_log10pi)
524 jmc 1.6 C The following form does the same, but is faster
525 jmc 1.21 mm_pi = EXP(mm_log10pi*lnTEN)
526     C over ice specific humidity
527 mlosch 1.5 qhice(I,J) = bb1*mm_pi/(Ppascals- (1.0 _d 0 - bb1) * mm_pi)
528     F_lh(I,J) = D1I * UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
529 ifenty 1.16 #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
530 jmc 1.21 IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
531 ifenty 1.16 F_lh(I,J) = F_lh_max(I,J)
532 jmc 1.21 ENDIF
533     #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
534 ifenty 1.16
535 jmc 1.22 F_c(I,J) = effConduct(I,J) * (TB(I,J) - t1)
536 mlosch 1.18 F_lwu(I,J) = t4 * D3(I,J)
537 mlosch 1.5 F_sens(I,J) = D1 * UG(I,J) * (t1 - atempLoc(I,J))
538 jmc 1.1
539 jmc 1.21 C The flux between the ice/snow surface and the atmosphere.
540 jmc 1.3 C (excludes upward conductive fluxes)
541 jmc 1.21 F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
542     & + F_sens(I,J) + F_lh(I,J)
543 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
544 jmc 1.1
545 gforget 1.9 #ifdef SEAICE_MODIFY_GROWTH_ADJ
546     Cgf no additional dependency through solver, snow, etc.
547 jmc 1.21 IF ( SEAICEadjMODE.GE.2 ) THEN
548 gforget 1.9 CALL ZERO_ADJ_1D( 1, TSURF(I,J,bi,bj), myThid)
549     t1 = TSURF(I,J,bi,bj)
550     t2 = t1*t1
551     t3 = t2*t1
552     t4 = t2*t2
553     qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
554    
555 mlosch 1.18 A1(I,J)=0.3 _d 0 *SWDOWN(I,J,bi,bj)+lwdownLoc(I,J)
556 gforget 1.9 & +D1*UG(I,J)*atempLoc(I,J)+D1I*UG(I,J)*AQH(I,J,bi,bj)
557 jmc 1.21 A2(I,J)= D1*UG(I,J)*t1+D1I*UG(I,J)*qhice(I,J)+D3(I,J)*t4
558 gforget 1.9
559 jmc 1.21 F_ia(I,J)=-A1(I,J)+A2(I,J)
560     IcePenetSW(I,J)= 0. _d 0
561     ENDIF
562     #endif /* SEAICE_MODIFY_GROWTH_ADJ */
563    
564     C Fresh water flux (kg/m^2/s) from latent heat of sublimation.
565     C F_lh is positive upward (sea ice looses heat) and FWsublim
566     C is also positive upward (atmosphere gains freshwater)
567     FWsublim(I,J) = F_lh(I,J)/lhSublim
568 gforget 1.9
569 jmc 1.21 #ifdef SEAICE_DEBUG
570 mlosch 1.5 C Caclulate the net ice-ocean and ice-atmosphere fluxes
571 jmc 1.22 IF (F_c(I,J) .GT. 0.0 _d 0) THEN
572     F_io_net(I,J) = F_c(I,J)
573 mlosch 1.5 F_ia_net(I,J) = 0.0 _d 0
574     ELSE
575     F_io_net(I,J) = 0.0 _d 0
576     F_ia_net(I,J) = F_ia(I,J)
577     ENDIF !/* conductive fluxes up or down */
578 jmc 1.1
579 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
580 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
581 mlosch 1.5 print '(A)','----------------------------------------'
582     print '(A,i6)','ibi complete ', myIter
583     print '(A,4(1x,D24.15))',
584     & 'ibi T(SURF, surfLoc,atmos) ',
585     & TSURF(I,J,bi,bj), tsurfLoc(I,J),atempLoc(I,J)
586     print '(A,4(1x,D24.15))',
587     & 'ibi LWL ', lwdownLoc(I,J)
588     print '(A,4(1x,D24.15))',
589     & 'ibi QSW(Total, Penetrating)',
590 jmc 1.21 & SWDOWN(I,J,bi,bj), IcePenetSW(I,J)
591 mlosch 1.5 print '(A,4(1x,D24.15))',
592     & 'ibi qh(ATM ICE) ',
593     & AQH(I,J,bi,bj),qhice(I,J)
594 ifenty 1.16 #ifndef SEAICE_SOLVE4TEMP_LEGACY
595     print '(A,4(1x,D24.15))',
596     & 'ibi F(lwd,swi,lwu) ',
597 jmc 1.21 & -lwdownLoc(I,J), -absorbedSW(I,J), F_lwu(I,J)
598 ifenty 1.16 print '(A,4(1x,D24.15))',
599     & 'ibi F(c,lh,sens) ',
600     & F_c(I,J), F_lh(I,J), F_sens(I,J)
601     #ifdef SEAICE_ADD_SUBLIMATION_TO_FWBUDGET
602     IF (F_lh_max(I,J) .GT. ZERO) THEN
603     print '(A,4(1x,D24.15))',
604     & 'ibi F_lh_max, F_lh/lhmax) ',
605     & F_lh_max(I,J), F_lh(I,J)/ F_lh_max(I,J)
606 jmc 1.19 ELSE
607 ifenty 1.16 print '(A,4(1x,D24.15))',
608     & 'ibi F_lh_max = ZERO! '
609     ENDIF
610     print '(A,4(1x,D24.15))',
611     & 'ibi FWsub, FWsubm*dT/rhoI ',
612     & FWsublim(I,J),
613     & FWsublim(I,J)*SEAICE_deltaTtherm/SEAICE_rhoICE
614 jmc 1.21 #endif /* SEAICE_ADD_SUBLIMATION_TO_FWBUDGET */
615     #endif /* SEAICE_SOLVE4TEMP_LEGACY */
616 mlosch 1.5 print '(A,4(1x,D24.15))',
617     & 'ibi F_ia, F_ia_net, F_c ',
618 jmc 1.6 #ifdef SEAICE_SOLVE4TEMP_LEGACY
619 jmc 1.22 & -A1(I,J)+A2(I,J), -A1(I,J)+A2(I,J)-F_c(I,J), F_c(I,J)
620 jmc 1.21 #else /* SEAICE_SOLVE4TEMP_LEGACY */
621     & F_ia(I,J), F_ia_net(I,J), F_c(I,J)
622 jmc 1.6 #endif /* SEAICE_SOLVE4TEMP_LEGACY */
623 mlosch 1.5 print '(A)','----------------------------------------'
624     ENDIF
625 jmc 1.2 #endif /* SEAICE_DEBUG */
626 jmc 1.6
627 mlosch 1.5 ENDIF !/* iceOrNot */
628     ENDDO !/* i */
629 jmc 1.1 ENDDO !/* j */
630    
631 jmc 1.19 #endif /* ALLOW_ATM_TEMP && ALLOW_DOWNWARD_RADIATION */
632     RETURN
633 jmc 1.1 END
  ViewVC Help
Powered by ViewVC 1.1.22