/[MITgcm]/MITgcm/pkg/thsice/thsice_step_fwd.F
ViewVC logotype

Contents of /MITgcm/pkg/thsice/thsice_step_fwd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.22 - (show annotations) (download)
Thu Apr 19 04:31:28 2007 UTC (17 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59
Changes since 1.21: +2 -1 lines
fix albedo & snow-precip diagnostics (when using advection).

1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.21 2007/04/17 23:42:33 heimbach Exp $
2 C $Name: $
3
4 #include "THSICE_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: THSICE_STEP_FWD
8 C !INTERFACE:
9 SUBROUTINE THSICE_STEP_FWD(
10 I bi, bj, iMin, iMax, jMin, jMax,
11 I prcAtm,
12 I myTime, myIter, myThid )
13 C !DESCRIPTION: \bv
14 C *==========================================================*
15 C | S/R THSICE_STEP_FWD
16 C | o Step Forward Therm-SeaIce model.
17 C *==========================================================*
18 C \ev
19
20 C !USES:
21 IMPLICIT NONE
22
23 C === Global variables ===
24 #include "SIZE.h"
25 #include "EEPARAMS.h"
26 #include "PARAMS.h"
27 #include "FFIELDS.h"
28 #include "THSICE_SIZE.h"
29 #include "THSICE_PARAMS.h"
30 #include "THSICE_VARS.h"
31 #include "THSICE_TAVE.h"
32 #include "THSICE_2DYN.h"
33 #ifdef ALLOW_AUTODIFF_TAMC
34 # include "tamc.h"
35 # include "tamc_keys.h"
36 #endif
37
38 INTEGER siLo, siHi, sjLo, sjHi
39 PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
40 PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
41
42 C !INPUT/OUTPUT PARAMETERS:
43 C === Routine arguments ===
44 C- input:
45 C bi,bj :: tile indices
46 C iMin,iMax :: computation domain: 1rst index range
47 C jMin,jMax :: computation domain: 2nd index range
48 C prcAtm :: total precip from the atmosphere [kg/m2/s]
49 C myTime :: current Time of simulation [s]
50 C myIter :: current Iteration number in simulation
51 C myThid :: my Thread Id number
52 C-- Use fluxes hold in commom blocks
53 C- input:
54 C icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
55 C icFlxAtm :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
56 C icFrwAtm :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
57 C- output
58 C icFlxAtm :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
59 C icFrwAtm :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
60 INTEGER bi,bj
61 INTEGER iMin, iMax
62 INTEGER jMin, jMax
63 _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 _RL myTime
65 INTEGER myIter
66 INTEGER myThid
67 CEOP
68
69 #ifdef ALLOW_THSICE
70 C !LOCAL VARIABLES:
71 C === Local variables ===
72 C iceFrac :: fraction of grid area covered in ice
73 C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2]
74 C frw2oc :: fresh-water flux from the ice to the ocean
75 C fsalt :: mass salt flux to the ocean
76 C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
77 C tFrzOce :: sea-water freezing temperature [oC] (function of S)
78 C isIceFree :: true for ice-free grid-cell that remains ice-free
79 C ageFac :: snow aging factor [1]
80 C snowFac :: snowing refreshing-age factor [units of 1/snowPr]
81 LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82 _RL iceFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83 _RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84 _RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85 _RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86 _RL tFrzOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 _RL frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL ageFac
89 _RL snowFac
90 _RL cphm
91 _RL opFrac, icFrac
92 #ifdef ALLOW_DIAGNOSTICS
93 _RL tmpFac
94 #endif
95 INTEGER i,j
96 LOGICAL dBugFlag
97
98 C- define grid-point location where to print debugging values
99 #include "THSICE_DEBUG.h"
100
101 1010 FORMAT(A,1P4E14.6)
102
103 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
104
105 #ifdef ALLOW_AUTODIFF_TAMC
106 act1 = bi - myBxLo(myThid)
107 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
108 act2 = bj - myByLo(myThid)
109 max2 = myByHi(myThid) - myByLo(myThid) + 1
110 act3 = myThid - 1
111 max3 = nTx*nTy
112 act4 = ikey_dynamics - 1
113 iicekey = (act1 + 1) + act2*max1
114 & + act3*max1*max2
115 & + act4*max1*max2*max3
116 #endif /* ALLOW_AUTODIFF_TAMC */
117
118 C- Initialise
119 dBugFlag = debugLevel.GE.debLevB
120 DO j = 1-OLy, sNy+OLy
121 DO i = 1-OLx, sNx+OLx
122 isIceFree(i,j) = .FALSE.
123 saltFlux(i,j,bi,bj) = 0. _d 0
124 #ifdef ALLOW_AUTODIFF_TAMC
125 iceFrac(i,j) = 0.
126 #endif
127 ENDDO
128 ENDDO
129
130 ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
131 snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
132
133 #ifdef ALLOW_AUTODIFF_TAMC
134 CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
135 #endif
136 DO j = jMin, jMax
137 DO i = iMin, iMax
138 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
139 C-- Snow aging :
140 snowAge(i,j,bi,bj) = thSIce_deltaT
141 & + snowAge(i,j,bi,bj)*ageFac
142 IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
143 & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
144 & * EXP( - snowFac*snowPrc(i,j,bi,bj) )
145 c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
146 c & /hNewSnowAge )
147 C-------
148 C note: Any flux of mass (here fresh water) that enter or leave the system
149 C with a non zero energy HAS TO be counted: add snow precip.
150 icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
151 & - Lfresh*snowPrc(i,j,bi,bj)
152 C--
153 ENDIF
154 ENDDO
155 ENDDO
156
157 #ifdef ALLOW_DIAGNOSTICS
158 IF ( useDiagnostics ) THEN
159 tmpFac = 1. _d 0
160 CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
161 CALL DIAGNOSTICS_FRACT_FILL(
162 I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
163 I 0,1,1,bi,bj,myThid)
164 CALL DIAGNOSTICS_FRACT_FILL(
165 I siceAlb, iceMask,tmpFac,1,'SIalbedo',
166 I 0,1,1,bi,bj,myThid)
167 ENDIF
168 #endif /* ALLOW_DIAGNOSTICS */
169 DO j = jMin, jMax
170 DO i = iMin, iMax
171 siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
172 ENDDO
173 ENDDO
174
175 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
176 C part.2 : ice-covered fraction ;
177 C change in ice/snow thickness and ice-fraction
178 C note: can only reduce the ice-fraction but not increase it.
179 C-------
180 DO j = jMin, jMax
181 DO i = iMin, iMax
182
183 tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
184 cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
185 frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
186 & * cphm/ocean_deltaT
187 iceFrac(i,j) = iceMask(i,j,bi,bj)
188 flx2oc(i,j) = icFlxSW(i,j,bi,bj)
189 C-------
190 #ifdef ALLOW_DBUG_THSICE
191 IF ( dBug(i,j,bi,bj) ) THEN
192 IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
193 WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
194 WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
195 & iceFrac(i,j), iceHeight(i,j,bi,bj),
196 & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
197 WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
198 & tOceMxL(i,j,bi,bj), tFrzOce(i,j),
199 & frzmltMxL(i,j), Qnet(i,j,bi,bj)
200 ENDIF
201 IF (iceFrac(i,j).GT.0.)
202 & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
203 & iceFrac(i,j), icFlxAtm(i,j,bi,bj),
204 & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
205 ENDIF
206 #endif
207 ENDDO
208 ENDDO
209
210 #ifdef ALLOW_AUTODIFF_TAMC
211 CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
212 #endif
213
214 CALL THSICE_CALC_THICKN(
215 I bi, bj, siLo, siHi, sjLo, sjHi,
216 I iMin,iMax, jMin,jMax, dBugFlag,
217 I iceMask(siLo,sjLo,bi,bj), tFrzOce,
218 I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
219 I snowPrc(siLo,sjLo,bi,bj), prcAtm,
220 I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
221 U iceFrac, iceHeight(siLo,sjLo,bi,bj),
222 U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
223 U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
224 U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
225 O frw2oc, fsalt,
226 I myTime, myIter, myThid )
227
228 C-- Net fluxes :
229 DO j = jMin, jMax
230 DO i = iMin, iMax
231 #ifdef ALLOW_AUTODIFF_TAMC
232 ikey_1 = i
233 & + sNx*(j-1)
234 & + sNx*sNy*act1
235 & + sNx*sNy*max1*act2
236 & + sNx*sNy*max1*max2*act3
237 & + sNx*sNy*max1*max2*max3*act4
238 #endif /* ALLOW_AUTODIFF_TAMC */
239 C--
240 #ifdef ALLOW_AUTODIFF_TAMC
241 CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
242 #endif
243 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
244 C- weighted average net fluxes:
245 #ifdef ALLOW_AUTODIFF_TAMC
246 CADJ STORE fsalt(i,j) = comlev1_thsice_1, key=ikey_1
247 CADJ STORE flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
248 CADJ STORE frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
249 CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
250 #endif
251 icFrac = iceMask(i,j,bi,bj)
252 opFrac= 1. _d 0-icFrac
253 icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
254 & - opFrac*Qnet(i,j,bi,bj)
255 icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
256 & + opFrac*rhofw*EmPmR(i,j,bi,bj)
257 Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
258 EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
259 & + opFrac*EmPmR(i,j,bi,bj)
260 saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
261
262 #ifdef ALLOW_DBUG_THSICE
263 IF (dBug(i,j,bi,bj)) WRITE(6,1010)
264 & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
265 & iceFrac(i,j), iceHeight(i,j,bi,bj),
266 & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
267 #endif
268
269 ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
270 icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
271 icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
272 ELSE
273 icFlxAtm(i,j,bi,bj) = 0. _d 0
274 icFrwAtm(i,j,bi,bj) = 0. _d 0
275 ENDIF
276 ENDDO
277 ENDDO
278
279 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
280 C part.3 : freezing of sea-water
281 C over ice-free fraction and what is left from ice-covered fraction
282 C-------
283 CALL THSICE_EXTEND(
284 I bi, bj, siLo, siHi, sjLo, sjHi,
285 I iMin,iMax, jMin,jMax, dBugFlag,
286 I frzmltMxL, tFrzOce,
287 I tOceMxL(siLo,sjLo,bi,bj),
288 U iceFrac, iceHeight(siLo,sjLo,bi,bj),
289 U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
290 U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
291 U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
292 O flx2oc, frw2oc, fsalt,
293 I myTime, myIter, myThid )
294
295 #ifdef ALLOW_AUTODIFF_TAMC
296 CADJ STORE snowHeight(:,:,bi,bj) =
297 CADJ & comlev1_bibj, key=iicekey, byte=isbyte
298 #endif
299 DO j = jMin, jMax
300 DO i = iMin, iMax
301 IF (frzmltMxL(i,j).GT.0. _d 0) THEN
302 C-- Net fluxes :
303 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
304 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
305 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
306
307 #ifdef ALLOW_DBUG_THSICE
308 IF (dBug(i,j,bi,bj)) WRITE(6,1010)
309 & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
310 & iceFrac(i,j), iceHeight(i,j,bi,bj),
311 & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
312 #endif
313 ENDIF
314
315 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
316 & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
317 & .AND. iceFrac(i,j) .LE.0. _d 0
318 IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
319 iceMask(i,j,bi,bj)=iceFrac(i,j)
320 IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
321 & snowAge(i,j,bi,bj) = 0. _d 0
322 ELSE
323 iceMask(i,j,bi,bj) = 0. _d 0
324 iceHeight(i,j,bi,bj)= 0. _d 0
325 snowHeight(i,j,bi,bj)=0. _d 0
326 snowAge(i,j,bi,bj) = 0. _d 0
327 Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
328 Tice1(i,j,bi,bj) = 0. _d 0
329 Tice2(i,j,bi,bj) = 0. _d 0
330 Qice1(i,j,bi,bj) = Lfresh
331 Qice2(i,j,bi,bj) = Lfresh
332 ENDIF
333 ENDDO
334 ENDDO
335
336 #ifdef ATMOSPHERIC_LOADING
337 # ifdef ALLOW_AUTODIFF_TAMC
338 CADJ STORE snowHeight(:,:,bi,bj) =
339 CADJ & comlev1_bibj, key=iicekey, byte=isbyte
340 # endif
341 DO j = jMin, jMax
342 DO i = iMin, iMax
343 C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
344 sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
345 & + iceHeight(i,j,bi,bj)*rhoi
346 & )*iceMask(i,j,bi,bj)
347 ENDDO
348 ENDDO
349 #endif
350
351 IF ( thSIceAdvScheme.GT.0 ) THEN
352 C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
353 DO j = jMin, jMax
354 DO i = iMin, iMax
355 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
356 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
357 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw
358 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
359 ENDIF
360 ENDDO
361 ENDDO
362 ENDIF
363
364 #ifdef ALLOW_BULK_FORCE
365 IF ( useBulkForce ) THEN
366 CALL BULKF_FLUX_ADJUST(
367 I bi, bj, iMin, iMax, jMin, jMax,
368 I isIceFree, myTime, myIter, myThid )
369 ENDIF
370 #endif /* ALLOW_BULK_FORCE */
371
372 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
373 #endif /* ALLOW_THSICE */
374
375 RETURN
376 END

  ViewVC Help
Powered by ViewVC 1.1.22