/[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.30 - (hide annotations) (download)
Sat Feb 11 03:35:01 2012 UTC (12 years, 3 months ago) by gforget
Branch: MAIN
Changes since 1.29: +25 -35 lines 
- unified freezing point treatment : the old SEAICE_VARIABLE_FREEZING_POINT
  is now the default, and the old default constant freezing point is recovered with
  SEAICE_tempFrz0    = -1.96,
  SEAICE_dTempFrz_dS = 0.,
- retire SEAICE_freeze that was the old way of specifycing the constant freezing point.
- remove ALLOW_SEAICE_FLOODING brackets; run time switch is already there.
- bug fix (thanks to M. Losch) : the sublimation term that was missing is now
  activated. To reproduce old results that had this bug define SEAICE_DISABLE_SUBLIM.
- bug fix (silly me) : flooding sign term (seaice_growth.F r 1.149line 1211)
  was wrong for SIsal0.NE.0. Changes cs32x15 results (see upcoming checkin for detail).
- bug fix (thanks to O. Jahn) : area loss for melting  got messed up in
  seaice_growth r149 for legacy branch (no results change).
- introduce SEAICE_CAP_SUBLIM : caps sublimation heat flux in solve4temp (code from I. Fenty).
- results did not change except for global_ocean.cs32x15, mostly due to the switch
  to variable freezing point (see upcoming checkin of results for details).
1 gforget 1.30 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_solve4temp.F,v 1.29 2012/02/07 14:01:10 mlosch 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 gforget 1.30 #ifdef SEAICE_CAP_SUBLIM
15 ifenty 1.16 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     #include "DYNVARS.h"
43     #ifdef ALLOW_EXF
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 jmc 1.21 C !INPUT PARAMETERS:
51     C UG :: atmospheric wind speed (m/s)
52 jmc 1.1 C HICE_ACTUAL :: actual ice thickness
53     C HSNOW_ACTUAL :: actual snow thickness
54 jmc 1.21 C TSURF :: surface temperature of ice/snow in Kelvin
55     C bi,bj :: tile indices
56     C myTime :: current time in simulation
57     C myIter :: iteration number in simulation
58     C myThid :: my Thread Id number
59     C !OUTPUT PARAMETERS:
60     C TSURF :: updated surface temperature of ice/snow in Kelvin
61     C F_ia :: upward seaice/snow surface heat flux to atmosphere (W/m^2)
62     C IcePenetSW :: short wave heat flux transmitted through ice (+=upward)
63     C FWsublim :: fresh water (mass) flux due to sublimation (+=up)(kg/m^2/s)
64 jmc 1.24 C---- Notes:
65     C 1) should add IcePenetSW to F_ia to get the net surface heat flux
66     C from the atmosphere (IcePenetSW not currently included in F_ia)
67     C 2) since zero ice/snow heat capacity is assumed, all the absorbed Short
68     C -Wave is used to warm the ice/snow surface (heating profile ignored).
69     C----------
70 jmc 1.21 _RL UG (1:sNx,1:sNy)
71     _RL HICE_ACTUAL (1:sNx,1:sNy)
72     _RL HSNOW_ACTUAL(1:sNx,1:sNy)
73 gforget 1.30 #ifdef SEAICE_CAP_SUBLIM
74 jmc 1.21 _RL F_lh_max (1:sNx,1:sNy)
75 ifenty 1.16 #endif
76 jmc 1.21 _RL TSURF (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
77     _RL F_ia (1:sNx,1:sNy)
78     _RL IcePenetSW (1:sNx,1:sNy)
79     _RL FWsublim (1:sNx,1:sNy)
80 jmc 1.1 INTEGER bi, bj
81     _RL myTime
82     INTEGER myIter, myThid
83 jmc 1.19 CEOP
84 jmc 1.1
85 jmc 1.19 #if defined(ALLOW_ATM_TEMP) && defined(ALLOW_DOWNWARD_RADIATION)
86 jmc 1.1 C !LOCAL VARIABLES:
87     C === Local variables ===
88 jmc 1.3 C i, j :: Loop counters
89 gforget 1.30 C kSurface :: vertical index of surface layer
90 jmc 1.1 INTEGER i, j
91 gforget 1.30 INTEGER kSurface
92 jmc 1.1 INTEGER ITER
93 gforget 1.30 C tempFrz :: ocean temperature in contact with ice (=seawater freezing point) (K)
94     _RL tempFrz (1:sNx,1:sNy)
95 mlosch 1.18 _RL D1, D1I
96     _RL D3(1:sNx,1:sNy)
97 mlosch 1.5 _RL TMELT, XKI, XKS, HCUT, XIO
98 jmc 1.27 C SurfMeltTemp :: Temp (K) above which wet-albedo values are used
99 mlosch 1.5 _RL SurfMeltTemp
100 jmc 1.21 C effConduct :: effective conductivity of combined ice and snow
101 mlosch 1.5 _RL effConduct(1:sNx,1:sNy)
102 jmc 1.21 C lhSublim :: latent heat of sublimation (SEAICE_lhEvap + SEAICE_lhFusion)
103     _RL lhSublim
104     C t1,t2,t3,t4 :: powers of temperature
105     _RL t1, t2, t3, t4
106 jmc 1.1
107 jmc 1.24 C- Constants to calculate Saturation Vapor Pressure
108 jmc 1.26 C Maykut Polynomial Coeff. for Sat. Vapor Press
109 jmc 1.21 _RL C1, C2, C3, C4, C5, QS1
110 jmc 1.24 C Extended temp-range expon. relation Coeff. for Sat. Vapor Press
111 jmc 1.21 _RL lnTEN
112 jmc 1.1 _RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
113     C specific humidity at ice surface variables
114 jmc 1.21 _RL mm_pi,mm_log10pi
115 jmc 1.1
116 jmc 1.22 C F_c :: conductive heat flux through seaice+snow (+=upward)
117 jmc 1.26 C F_lwu :: upward long-wave surface heat flux (+=upward)
118     C F_sens :: sensible surface heat flux (+=upward)
119 jmc 1.21 C F_lh :: latent heat flux (sublimation) (+=upward)
120 jmc 1.24 C qhice :: saturation vapor pressure of snow/ice surface
121     C dqh_dTs :: derivative of qhice w.r.t snow/ice surf. temp
122     C dFia_dTs :: derivative of surf heat flux (F_ia) w.r.t surf. temp
123 jmc 1.26 _RL F_c (1:sNx,1:sNy)
124 jmc 1.21 _RL F_lwu (1:sNx,1:sNy)
125     _RL F_sens (1:sNx,1:sNy)
126     _RL F_lh (1:sNx,1:sNy)
127     _RL qhice (1:sNx,1:sNy)
128 jmc 1.24 _RL dqh_dTs (1:sNx,1:sNy)
129 jmc 1.26 _RL dFia_dTs (1:sNx,1:sNy)
130 jmc 1.21 _RL absorbedSW (1:sNx,1:sNy)
131     _RL penetSWFrac
132 jmc 1.25 _RL delTsurf
133 jmc 1.21
134     C local copies of global variables
135     _RL tsurfLoc (1:sNx,1:sNy)
136 jmc 1.25 _RL tsurfPrev (1:sNx,1:sNy)
137 jmc 1.21 _RL atempLoc (1:sNx,1:sNy)
138     _RL lwdownLoc (1:sNx,1:sNy)
139     _RL ALB (1:sNx,1:sNy)
140     _RL ALB_ICE (1:sNx,1:sNy)
141     _RL ALB_SNOW (1:sNx,1:sNy)
142     C iceOrNot :: this is HICE_ACTUAL.GT.0.
143     LOGICAL iceOrNot(1:sNx,1:sNy)
144     #ifdef SEAICE_DEBUG
145     C F_io_net :: upward conductive heat flux through seaice+snow
146     C F_ia_net :: net heat flux divergence at the sea ice/snow surface:
147     C includes ice conductive fluxes and atmospheric fluxes (W/m^2)
148 jmc 1.25 _RL F_io_net
149     _RL F_ia_net
150 jmc 1.21 #endif /* SEAICE_DEBUG */
151 ifenty 1.14
152 jmc 1.21 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
153 jmc 1.1
154 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
155     CADJ INIT comlev1_solve4temp = COMMON, sNx*sNy*NMAX_TICE
156     #endif /* ALLOW_AUTODIFF_TAMC */
157    
158 jmc 1.24 C- MAYKUT 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 QS1=0.622 _d +00/1013.0 _d +00
165 jmc 1.24 C- Extended temp-range expon. relation Coeff. for Sat. Vapor Press
166 jmc 1.21 lnTEN = LOG(10.0 _d 0)
167 jmc 1.1 aa1 = 2663.5 _d 0
168     aa2 = 12.537 _d 0
169     bb1 = 0.622 _d 0
170 mlosch 1.5 bb2 = 1.0 _d 0 - bb1
171 jmc 1.1 Ppascals = 100000. _d 0
172 mlosch 1.5 C cc0 = TEN ** aa2
173 jmc 1.21 cc0 = EXP(aa2*lnTEN)
174 mlosch 1.5 cc1 = cc0*aa1*bb1*Ppascals*lnTEN
175 jmc 1.1 cc2 = cc0*bb2
176    
177 gforget 1.30 IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
178     kSurface = Nr
179     ELSE
180     kSurface = 1
181     ENDIF
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     C MELTING TEMPERATURE OF ICE
191 mlosch 1.5 TMELT = celsius2K
192 jmc 1.1
193     C ICE CONDUCTIVITY
194     XKI=SEAICE_iceConduct
195    
196     C SNOW CONDUCTIVITY
197     XKS=SEAICE_snowConduct
198    
199     C CUTOFF SNOW THICKNESS
200 jmc 1.21 C Snow-Thickness above HCUT: SW optically thick snow (=> snow-albedo).
201     C Snow-Thickness below HCUT: linear transition to ice-albedo
202 jmc 1.24 HCUT = SEAICE_snowThick
203    
204     C PENETRATION SHORTWAVE RADIATION FACTOR
205     XIO=SEAICE_shortwave
206    
207 jmc 1.27 C Temperature Threshold for wet-albedo:
208     SurfMeltTemp = TMELT + SEAICE_wetAlbTemp
209 jmc 1.24 C old SOLVE4TEMP_LEGACY setting, consistent with former celsius2K value:
210     c TMELT = 273.16 _d +00
211     c SurfMeltTemp = 273.159 _d +00
212 jmc 1.1
213 jmc 1.3 C Initialize variables
214 jmc 1.1 DO J=1,sNy
215 mlosch 1.5 DO I=1,sNx
216     C HICE_ACTUAL is modified in this routine, but at the same time
217     C used to decided where there is ice, therefore we save this information
218     C here in a separate array
219 jmc 1.21 iceOrNot (I,J) = HICE_ACTUAL(I,J) .GT. 0. _d 0
220     IcePenetSW(I,J) = 0. _d 0
221     absorbedSW(I,J) = 0. _d 0
222 mlosch 1.5 qhice (I,J) = 0. _d 0
223 jmc 1.24 dqh_dTs (I,J) = 0. _d 0
224 mlosch 1.5 F_ia (I,J) = 0. _d 0
225 mlosch 1.10 F_lh (I,J) = 0. _d 0
226 mlosch 1.5 F_lwu (I,J) = 0. _d 0
227     F_sens (I,J) = 0. _d 0
228 jmc 1.26 C Make a local copy of LW, surface & atmospheric temperatures
229 mlosch 1.5 tsurfLoc (I,J) = TSURF(I,J,bi,bj)
230 jmc 1.26 c tsurfLoc (I,J) = MIN( celsius2K+MAX_TICE, TSURF(I,J,bi,bj) )
231     lwdownLoc(I,J) = MAX( MIN_LWDOWN, LWDOWN(I,J,bi,bj) )
232 jmc 1.24 atempLoc (I,J) = MAX( celsius2K+MIN_ATEMP, ATEMP(I,J,bi,bj) )
233 jmc 1.1
234 gforget 1.30 c FREEZING TEMP. OF SEA WATER (K)
235     tempFrz(I,J) = SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
236     & + SEAICE_tempFrz0 + celsius2K
237    
238 jmc 1.26 C Now determine fixed (relative to tsurf) forcing term in heat budget
239    
240 mlosch 1.18 IF(HSNOW_ACTUAL(I,J).GT.0.0) THEN
241 jmc 1.21 C Stefan-Boltzmann constant times emissivity
242 mlosch 1.18 D3(I,J)=SEAICE_snow_emiss*SEAICE_boltzmann
243     #ifdef EXF_LWDOWN_WITH_EMISSIVITY
244     C This is now [(1-emiss)*lwdown - lwdown]
245 jmc 1.26 lwdownLoc(I,J) = SEAICE_snow_emiss*lwdownLoc(I,J)
246 mlosch 1.18 #else /* use the old hard wired inconsistent value */
247 jmc 1.26 lwdownLoc(I,J) = 0.97 _d 0*lwdownLoc(I,J)
248 mlosch 1.18 #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
249     ELSE
250 jmc 1.21 C Stefan-Boltzmann constant times emissivity
251 mlosch 1.18 D3(I,J)=SEAICE_ice_emiss*SEAICE_boltzmann
252     #ifdef EXF_LWDOWN_WITH_EMISSIVITY
253     C This is now [(1-emiss)*lwdown - lwdown]
254 jmc 1.26 lwdownLoc(I,J) = SEAICE_ice_emiss*lwdownLoc(I,J)
255 mlosch 1.18 #else /* use the old hard wired inconsistent value */
256 jmc 1.26 lwdownLoc(I,J) = 0.97 _d 0*lwdownLoc(I,J)
257 mlosch 1.18 #endif /* EXF_LWDOWN_WITH_EMISSIVITY */
258     ENDIF
259 mlosch 1.5 ENDDO
260 jmc 1.1 ENDDO
261    
262     DO J=1,sNy
263 mlosch 1.5 DO I=1,sNx
264 jmc 1.1
265     C DECIDE ON ALBEDO
266 mlosch 1.5 IF ( iceOrNot(I,J) ) THEN
267 jmc 1.6
268 mlosch 1.5 IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 ) THEN
269     IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
270     ALB_ICE (I,J) = SEAICE_wetIceAlb_south
271     ALB_SNOW(I,J) = SEAICE_wetSnowAlb_south
272     ELSE ! no surface melting
273     ALB_ICE (I,J) = SEAICE_dryIceAlb_south
274     ALB_SNOW(I,J) = SEAICE_drySnowAlb_south
275     ENDIF
276     ELSE !/ Northern Hemisphere
277     IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
278     ALB_ICE (I,J) = SEAICE_wetIceAlb
279     ALB_SNOW(I,J) = SEAICE_wetSnowAlb
280     ELSE ! no surface melting
281     ALB_ICE (I,J) = SEAICE_dryIceAlb
282     ALB_SNOW(I,J) = SEAICE_drySnowAlb
283     ENDIF
284     ENDIF !/ Albedo for snow and ice
285    
286 jmc 1.21 C If actual snow thickness exceeds the cutoff thickness, use snow albedo
287 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. HCUT) THEN
288     ALB(I,J) = ALB_SNOW(I,J)
289 jmc 1.21 ELSEIF ( HCUT.LE.ZERO ) THEN
290     ALB(I,J) = ALB_ICE(I,J)
291 mlosch 1.5 ELSE
292 jmc 1.21 C otherwise, use linear transition between ice and snow albedo
293     ALB(I,J) = MIN( ALB_ICE(I,J) + HSNOW_ACTUAL(I,J)/HCUT
294     & *(ALB_SNOW(I,J) -ALB_ICE(I,J))
295     & , ALB_SNOW(I,J) )
296 mlosch 1.5 ENDIF
297    
298 jmc 1.21 C Determine the fraction of shortwave radiative flux remaining
299     C at ocean interface after scattering through the snow and ice.
300     C If snow is present, no radiation penetrates through snow+ice
301 mlosch 1.5 IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
302 jmc 1.21 penetSWFrac = 0.0 _d 0
303 mlosch 1.5 ELSE
304 jmc 1.21 penetSWFrac = XIO*EXP(-1.5 _d 0 * HICE_ACTUAL(I,J))
305 mlosch 1.5 ENDIF
306 jmc 1.21 C The shortwave radiative flux leaving ocean beneath ice (+=up).
307     IcePenetSW(I,J) = -(1.0 _d 0 - ALB(I,J))
308     & *penetSWFrac * SWDOWN(I,J,bi,bj)
309     C The shortwave radiative flux convergence in the seaice.
310     absorbedSW(I,J) = (1.0 _d 0 - ALB(I,J))
311     & *(1.0 _d 0 - penetSWFrac)* SWDOWN(I,J,bi,bj)
312 jmc 1.1
313 jmc 1.24 C The effective conductivity of the two-layer snow/ice system.
314 jmc 1.26 C Set a minimum sea ice thickness of 5 cm to bound
315 jmc 1.21 C the magnitude of conductive heat fluxes.
316     Cif * now taken care of by SEAICE_hice_reg in seaice_growth
317     c hice_tmp = max(HICE_ACTUAL(I,J),5. _d -2)
318 mlosch 1.5 effConduct(I,J) = XKI * XKS /
319 jmc 1.21 & (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
320 jmc 1.1
321     #ifdef SEAICE_DEBUG
322 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
323 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
324 mlosch 1.5 print '(A,i6)','-----------------------------------'
325     print '(A,i6)','ibi merged initialization ', myIter
326     print '(A,i6,4(1x,D24.15))',
327     & 'ibi iter, TSL, TS ',myIter,
328     & tsurfLoc(I,J), TSURF(I,J,bi,bj)
329     print '(A,i6,4(1x,D24.15))',
330     & 'ibi iter, TMELT ',myIter,TMELT
331     print '(A,i6,4(1x,D24.15))',
332     & 'ibi iter, HIA, EFKCON ',myIter,
333     & HICE_ACTUAL(I,J), effConduct(I,J)
334     print '(A,i6,4(1x,D24.15))',
335     & 'ibi iter, HSNOW ',myIter,
336     & HSNOW_ACTUAL(I,J), ALB(I,J)
337     print '(A,i6)','-----------------------------------'
338     print '(A,i6)','ibi energy balance iterat ', myIter
339     ENDIF
340 jmc 1.2 #endif /* SEAICE_DEBUG */
341 jmc 1.6
342 mlosch 1.5 ENDIF !/* iceOrNot */
343     ENDDO !/* i */
344     ENDDO !/* j */
345 jmc 1.21
346     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
347 mlosch 1.5 DO ITER=1,IMAX_TICE
348     DO J=1,sNy
349     DO I=1,sNx
350 mlosch 1.8 #ifdef ALLOW_AUTODIFF_TAMC
351     iicekey = I + sNx*(J-1) + (ITER-1)*sNx*sNy
352 jmc 1.26 CADJ STORE tsurfLoc(i,j) = comlev1_solve4temp,
353 mlosch 1.8 CADJ & key = iicekey, byte = isbyte
354     #endif /* ALLOW_AUTODIFF_TAMC */
355    
356 jmc 1.25 C- save tsurf from previous iter
357     tsurfPrev(I,J) = tsurfLoc(I,J)
358 mlosch 1.5 IF ( iceOrNot(I,J) ) THEN
359 jmc 1.1
360 mlosch 1.5 t1 = tsurfLoc(I,J)
361     t2 = t1*t1
362     t3 = t2*t1
363     t4 = t2*t2
364 jmc 1.1
365 jmc 1.24 C-- Calculate the specific humidity in the BL above the snow/ice
366 jmc 1.25 IF ( useMaykutSatVapPoly ) THEN
367 jmc 1.24 C- Use the Maykut polynomial
368 jmc 1.25 qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
369     dqh_dTs(I,J) = 0. _d 0
370     ELSE
371 jmc 1.24 C- Use exponential relation approx., more accurate at low temperatures
372 mlosch 1.5 C log 10 of the sat vap pressure
373 jmc 1.25 mm_log10pi = -aa1 / t1 + aa2
374 mlosch 1.5 C The saturation vapor pressure (SVP) in the surface
375     C boundary layer (BL) above the snow/ice.
376 jmc 1.25 c mm_pi = TEN **(mm_log10pi)
377 jmc 1.6 C The following form does the same, but is faster
378 jmc 1.25 mm_pi = EXP(mm_log10pi*lnTEN)
379     qhice(I,J) = bb1*mm_pi/( Ppascals -(1.0 _d 0 - bb1)*mm_pi )
380 jmc 1.21 C A constant for SVP derivative w.r.t TICE
381 jmc 1.25 c cc3t = TEN **(aa1 / t1)
382 jmc 1.21 C The following form does the same, but is faster
383 jmc 1.25 cc3t = EXP(aa1 / t1 * lnTEN)
384 jmc 1.21 C d(qh)/d(TICE)
385 jmc 1.25 dqh_dTs(I,J) = cc1*cc3t/((cc2-cc3t*Ppascals)**2 *t2)
386     ENDIF
387 jmc 1.1
388 mlosch 1.10 C Calculate the flux terms based on the updated tsurfLoc
389 gforget 1.30 F_c(I,J) = effConduct(I,J)*(tempFrz(I,J)-tsurfLoc(I,J))
390 mlosch 1.10 F_lh(I,J) = D1I*UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
391 gforget 1.30 #ifdef SEAICE_CAP_SUBLIM
392 jmc 1.21 C if the latent heat flux implied by tsurfLoc exceeds
393     C F_lh_max, cap F_lh and decouple the flux magnitude from TICE
394     IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
395 ifenty 1.16 F_lh(I,J) = F_lh_max(I,J)
396 jmc 1.24 dqh_dTs(I,J) = ZERO
397 jmc 1.21 ENDIF
398 gforget 1.30 #endif /* SEAICE_CAP_SUBLIM */
399 ifenty 1.16
400 jmc 1.21 F_lwu(I,J) = t4 * D3(I,J)
401 mlosch 1.5 F_sens(I,J)= D1 * UG(I,J) * (t1 - atempLoc(I,J))
402 jmc 1.21 F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
403 jmc 1.24 & + F_sens(I,J) + F_lh(I,J)
404 jmc 1.26 C d(F_ia)/d(Tsurf)
405     dFia_dTs(I,J) = 4.0 _d 0*D3(I,J)*t3 + D1*UG(I,J)
406     & + D1I*UG(I,J)*dqh_dTs(I,J)
407 jmc 1.1
408     #ifdef SEAICE_DEBUG
409 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
410 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
411 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
412     & 'ice-iter qhICE, ', ITER,qhIce(I,J)
413     print '(A,i6,4(1x,D24.15))',
414 jmc 1.24 & 'ice-iter dFiDTs1 F_ia ', ITER,
415     & dFia_dTs(I,J)+effConduct(I,J), F_ia(I,J)-F_c(I,J)
416 mlosch 1.5 ENDIF
417 jmc 1.2 #endif /* SEAICE_DEBUG */
418 jmc 1.1
419 jmc 1.26 C- Update tsurf as solution of : Fc = Fia + d/dT(Fia - Fc) *delta.tsurf
420 jmc 1.24 tsurfLoc(I,J) = tsurfLoc(I,J)
421     & + ( F_c(I,J)-F_ia(I,J) ) / ( effConduct(I,J)+dFia_dTs(I,J) )
422 jmc 1.1
423 jmc 1.26 IF ( useMaykutSatVapPoly ) THEN
424     tsurfLoc(I,J) = MAX( celsius2K+MIN_TICE, tsurfLoc(I,J) )
425     ENDIF
426 jmc 1.21 C If the search leads to tsurfLoc < 50 Kelvin, restart the search
427     C at tsurfLoc = TMELT. Note that one solution to the energy balance problem
428     C is an extremely low temperature - a temperature far below realistic values.
429 jmc 1.26 c IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) tsurfLoc(I,J) = TMELT
430     C Comments & code above not relevant anymore (from older version, when
431     C trying Maykut-Polynomial & dqh_dTs > 0 ?): commented out
432 jmc 1.21 tsurfLoc(I,J) = MIN( tsurfLoc(I,J), TMELT )
433 jmc 1.1
434     #ifdef SEAICE_DEBUG
435 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
436 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
437 mlosch 1.5 print '(A,i6,4(1x,D24.15))',
438     & 'ice-iter tsurfLc,|dif|', ITER,
439     & tsurfLoc(I,J),
440 jmc 1.21 & LOG10(ABS(tsurfLoc(I,J) - t1))
441 mlosch 1.5 ENDIF
442 jmc 1.2 #endif /* SEAICE_DEBUG */
443 jmc 1.1
444 mlosch 1.5 ENDIF !/* iceOrNot */
445     ENDDO !/* i */
446     ENDDO !/* j */
447     ENDDO !/* Iterations */
448 jmc 1.21 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
449    
450 mlosch 1.5 DO J=1,sNy
451     DO I=1,sNx
452     IF ( iceOrNot(I,J) ) THEN
453 jmc 1.1
454 jmc 1.21 C Save updated tsurf and finalize the flux terms
455     TSURF(I,J,bi,bj) = tsurfLoc(I,J)
456    
457 jmc 1.25 #ifdef SEAICE_MODIFY_GROWTH_ADJ
458     Cgf no additional dependency through solver, snow, etc.
459     IF ( SEAICEadjMODE.GE.2 ) THEN
460     CALL ZERO_ADJ_1D( 1, TSURF(I,J,bi,bj), myThid)
461 jmc 1.26 absorbedSW(I,J) = 0.3 _d 0 *SWDOWN(I,J,bi,bj)
462     IcePenetSW(I,J)= 0. _d 0
463     ENDIF
464     IF ( postSolvTempIter.EQ.2 .OR. SEAICEadjMODE.GE.2 ) THEN
465 jmc 1.25 t1 = TSURF(I,J,bi,bj)
466     #else /* SEAICE_MODIFY_GROWTH_ADJ */
467    
468     IF ( postSolvTempIter.EQ.2 ) THEN
469 mlosch 1.5 C Recalculate the fluxes based on the (possibly) adjusted TSURF
470 jmc 1.25 t1 = tsurfLoc(I,J)
471 jmc 1.26 #endif /* SEAICE_MODIFY_GROWTH_ADJ */
472 jmc 1.25 t2 = t1*t1
473     t3 = t2*t1
474     t4 = t2*t2
475 jmc 1.1
476 jmc 1.25 IF ( useMaykutSatVapPoly ) THEN
477     qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
478     ELSE
479 mlosch 1.5 C log 10 of the sat vap pressure
480 jmc 1.25 mm_log10pi = -aa1 / t1 + aa2
481 mlosch 1.5 C saturation vapor pressure
482 jmc 1.25 c mm_pi = TEN **(mm_log10pi)
483 jmc 1.6 C The following form does the same, but is faster
484 jmc 1.25 mm_pi = EXP(mm_log10pi*lnTEN)
485 jmc 1.21 C over ice specific humidity
486 jmc 1.25 qhice(I,J) = bb1*mm_pi/( Ppascals -(1.0 _d 0 - bb1)*mm_pi )
487     ENDIF
488 gforget 1.30 F_c(I,J) = effConduct(I,J) * (tempFrz(I,J) - t1)
489 jmc 1.25 F_lh(I,J) = D1I * UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
490 gforget 1.30 #ifdef SEAICE_CAP_SUBLIM
491 jmc 1.25 IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
492 ifenty 1.16 F_lh(I,J) = F_lh_max(I,J)
493 jmc 1.25 ENDIF
494 gforget 1.30 #endif /* SEAICE_CAP_SUBLIM */
495 jmc 1.25 F_lwu(I,J) = t4 * D3(I,J)
496     F_sens(I,J) = D1 * UG(I,J) * (t1 - atempLoc(I,J))
497 jmc 1.21 C The flux between the ice/snow surface and the atmosphere.
498 jmc 1.25 F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
499     & + F_sens(I,J) + F_lh(I,J)
500 jmc 1.1
501 jmc 1.25 ELSEIF ( postSolvTempIter.EQ.1 ) THEN
502     C Update fluxes (consistent with the linearized formulation)
503     delTsurf = tsurfLoc(I,J)-tsurfPrev(I,J)
504 gforget 1.30 F_c(I,J) = effConduct(I,J)*(tempFrz(I,J)-tsurfLoc(I,J))
505 jmc 1.25 F_ia(I,J) = F_ia(I,J) + dFia_dTs(I,J)*delTsurf
506     F_lh(I,J) = F_lh(I,J)
507     & + D1I*UG(I,J)*dqh_dTs(I,J)*delTsurf
508 jmc 1.26
509     c ELSEIF ( postSolvTempIter.EQ.0 ) THEN
510 jmc 1.25 C Take fluxes from last iteration
511 jmc 1.26
512 jmc 1.21 ENDIF
513    
514     C Fresh water flux (kg/m^2/s) from latent heat of sublimation.
515     C F_lh is positive upward (sea ice looses heat) and FWsublim
516     C is also positive upward (atmosphere gains freshwater)
517     FWsublim(I,J) = F_lh(I,J)/lhSublim
518 gforget 1.9
519 jmc 1.21 #ifdef SEAICE_DEBUG
520 jmc 1.26 C Calculate the net ice-ocean and ice-atmosphere fluxes
521 jmc 1.22 IF (F_c(I,J) .GT. 0.0 _d 0) THEN
522 jmc 1.25 F_io_net = F_c(I,J)
523     F_ia_net = 0.0 _d 0
524 mlosch 1.5 ELSE
525 jmc 1.25 F_io_net = 0.0 _d 0
526     F_ia_net = F_ia(I,J)
527 mlosch 1.5 ENDIF !/* conductive fluxes up or down */
528 jmc 1.1
529 jmc 1.21 IF ( (I .EQ. SEAICE_debugPointI) .AND.
530 ifenty 1.16 & (J .EQ. SEAICE_debugPointJ) ) THEN
531 mlosch 1.5 print '(A)','----------------------------------------'
532     print '(A,i6)','ibi complete ', myIter
533     print '(A,4(1x,D24.15))',
534     & 'ibi T(SURF, surfLoc,atmos) ',
535     & TSURF(I,J,bi,bj), tsurfLoc(I,J),atempLoc(I,J)
536     print '(A,4(1x,D24.15))',
537     & 'ibi LWL ', lwdownLoc(I,J)
538     print '(A,4(1x,D24.15))',
539     & 'ibi QSW(Total, Penetrating)',
540 jmc 1.21 & SWDOWN(I,J,bi,bj), IcePenetSW(I,J)
541 mlosch 1.5 print '(A,4(1x,D24.15))',
542     & 'ibi qh(ATM ICE) ',
543     & AQH(I,J,bi,bj),qhice(I,J)
544 ifenty 1.16 print '(A,4(1x,D24.15))',
545     & 'ibi F(lwd,swi,lwu) ',
546 jmc 1.21 & -lwdownLoc(I,J), -absorbedSW(I,J), F_lwu(I,J)
547 ifenty 1.16 print '(A,4(1x,D24.15))',
548     & 'ibi F(c,lh,sens) ',
549     & F_c(I,J), F_lh(I,J), F_sens(I,J)
550 gforget 1.30 #ifdef SEAICE_CAP_SUBLIM
551 ifenty 1.16 IF (F_lh_max(I,J) .GT. ZERO) THEN
552     print '(A,4(1x,D24.15))',
553     & 'ibi F_lh_max, F_lh/lhmax) ',
554     & F_lh_max(I,J), F_lh(I,J)/ F_lh_max(I,J)
555 jmc 1.19 ELSE
556 ifenty 1.16 print '(A,4(1x,D24.15))',
557     & 'ibi F_lh_max = ZERO! '
558     ENDIF
559     print '(A,4(1x,D24.15))',
560     & 'ibi FWsub, FWsubm*dT/rhoI ',
561     & FWsublim(I,J),
562     & FWsublim(I,J)*SEAICE_deltaTtherm/SEAICE_rhoICE
563 gforget 1.30 #endif /* SEAICE_CAP_SUBLIM */
564 mlosch 1.5 print '(A,4(1x,D24.15))',
565     & 'ibi F_ia, F_ia_net, F_c ',
566 jmc 1.25 & F_ia(I,J), F_ia_net, F_c(I,J)
567 mlosch 1.5 print '(A)','----------------------------------------'
568     ENDIF
569 jmc 1.2 #endif /* SEAICE_DEBUG */
570 jmc 1.6
571 mlosch 1.5 ENDIF !/* iceOrNot */
572     ENDDO !/* i */
573 jmc 1.1 ENDDO !/* j */
574    
575 jmc 1.19 #endif /* ALLOW_ATM_TEMP && ALLOW_DOWNWARD_RADIATION */
576     RETURN
577 jmc 1.1 END
  ViewVC Help
Powered by ViewVC 1.1.22