/[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.30 - (show annotations) (download)
Fri Dec 17 04:00:14 2010 UTC (13 years, 4 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62y, checkpoint62x
Changes since 1.29: +15 -15 lines
rename iicekey as ticekey to avoid conflict with pkg/seaice

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

  ViewVC Help
Powered by ViewVC 1.1.22