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