/[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.26 - (show annotations) (download)
Mon Oct 1 13:38:08 2007 UTC (16 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint60, checkpoint61, checkpoint62, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59k, checkpoint59j, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62l, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.25: +7 -7 lines
Change EmPmR units: from m/s to kg/m2/s

1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.25 2007/08/30 15:02:51 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 iicekey = (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=iicekey, byte=isbyte
145 #endif
146 DO j = jMin, jMax
147 DO i = iMin, iMax
148 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
149 C-- Snow aging :
150 snowAge(i,j,bi,bj) = thSIce_deltaT
151 & + snowAge(i,j,bi,bj)*ageFac
152 IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
153 & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
154 & * EXP( - snowFac*snowPrc(i,j,bi,bj) )
155 c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
156 c & /hNewSnowAge )
157 C-------
158 C note: Any flux of mass (here fresh water) that enter or leave the system
159 C with a non zero energy HAS TO be counted: add snow precip.
160 icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
161 & - Lfresh*snowPrc(i,j,bi,bj)
162 C--
163 ENDIF
164 ENDDO
165 ENDDO
166
167 #ifdef ALLOW_DIAGNOSTICS
168 IF ( useDiagnostics ) THEN
169 tmpFac = 1. _d 0
170 CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
171 CALL DIAGNOSTICS_FRACT_FILL(
172 I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
173 I 0,1,1,bi,bj,myThid)
174 CALL DIAGNOSTICS_FRACT_FILL(
175 I siceAlb, iceMask,tmpFac,1,'SIalbedo',
176 I 0,1,1,bi,bj,myThid)
177 ENDIF
178 #endif /* ALLOW_DIAGNOSTICS */
179 DO j = jMin, jMax
180 DO i = iMin, iMax
181 siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
182 ENDDO
183 ENDDO
184
185 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
186 C part.2 : ice-covered fraction ;
187 C change in ice/snow thickness and ice-fraction
188 C note: can only reduce the ice-fraction but not increase it.
189 C-------
190 DO j = jMin, jMax
191 DO i = iMin, iMax
192
193 tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
194 cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
195 frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
196 & * cphm/ocean_deltaT
197 iceFrac(i,j) = iceMask(i,j,bi,bj)
198 flx2oc(i,j) = icFlxSW(i,j,bi,bj)
199 C-------
200 #ifdef ALLOW_DBUG_THSICE
201 IF ( dBug(i,j,bi,bj) ) THEN
202 IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
203 WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
204 WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
205 & iceFrac(i,j), iceHeight(i,j,bi,bj),
206 & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
207 WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
208 & tOceMxL(i,j,bi,bj), tFrzOce(i,j),
209 & frzmltMxL(i,j), Qnet(i,j,bi,bj)
210 ENDIF
211 IF (iceFrac(i,j).GT.0.)
212 & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
213 & iceFrac(i,j), icFlxAtm(i,j,bi,bj),
214 & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
215 ENDIF
216 #endif
217 ENDDO
218 ENDDO
219
220 #ifdef ALLOW_AUTODIFF_TAMC
221 CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
222 #endif
223
224 CALL THSICE_CALC_THICKN(
225 I bi, bj, siLo, siHi, sjLo, sjHi,
226 I iMin,iMax, jMin,jMax, dBugFlag,
227 I iceMask(siLo,sjLo,bi,bj), tFrzOce,
228 I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
229 I snowPrc(siLo,sjLo,bi,bj), prcAtm,
230 I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
231 U iceFrac, iceHeight(siLo,sjLo,bi,bj),
232 U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
233 U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
234 U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
235 O frw2oc, fsalt,
236 I myTime, myIter, myThid )
237
238 C-- Net fluxes :
239 DO j = jMin, jMax
240 DO i = iMin, iMax
241 #ifdef ALLOW_AUTODIFF_TAMC
242 ikey_1 = i
243 & + sNx*(j-1)
244 & + sNx*sNy*act1
245 & + sNx*sNy*max1*act2
246 & + sNx*sNy*max1*max2*act3
247 & + sNx*sNy*max1*max2*max3*act4
248 #endif /* ALLOW_AUTODIFF_TAMC */
249 C--
250 #ifdef ALLOW_AUTODIFF_TAMC
251 CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
252 #endif
253 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
254 C- weighted average net fluxes:
255 #ifdef ALLOW_AUTODIFF_TAMC
256 CADJ STORE fsalt(i,j) = comlev1_thsice_1, key=ikey_1
257 CADJ STORE flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
258 CADJ STORE frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
259 CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
260 #endif
261 icFrac = iceMask(i,j,bi,bj)
262 opFrac= 1. _d 0-icFrac
263 #ifdef ALLOW_ATM2D
264 pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
265 pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
266 sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
267 #else
268 icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
269 & - opFrac*Qnet(i,j,bi,bj)
270 icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
271 & + opFrac*EmPmR(i,j,bi,bj)
272 Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
273 EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)
274 & + opFrac*EmPmR(i,j,bi,bj)
275 saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
276 #endif
277
278 #ifdef ALLOW_DBUG_THSICE
279 IF (dBug(i,j,bi,bj)) WRITE(6,1010)
280 & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
281 & iceFrac(i,j), iceHeight(i,j,bi,bj),
282 & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
283 #endif
284
285 ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN
286 icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
287 icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj)
288 ELSE
289 icFlxAtm(i,j,bi,bj) = 0. _d 0
290 icFrwAtm(i,j,bi,bj) = 0. _d 0
291 ENDIF
292 ENDDO
293 ENDDO
294
295 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
296 C part.3 : freezing of sea-water
297 C over ice-free fraction and what is left from ice-covered fraction
298 C-------
299 CALL THSICE_EXTEND(
300 I bi, bj, siLo, siHi, sjLo, sjHi,
301 I iMin,iMax, jMin,jMax, dBugFlag,
302 I frzmltMxL, tFrzOce,
303 I tOceMxL(siLo,sjLo,bi,bj),
304 U iceFrac, iceHeight(siLo,sjLo,bi,bj),
305 U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
306 U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
307 U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
308 O flx2oc, frw2oc, fsalt,
309 I myTime, myIter, myThid )
310
311 #ifdef ALLOW_AUTODIFF_TAMC
312 CADJ STORE snowHeight(:,:,bi,bj) =
313 CADJ & comlev1_bibj, key=iicekey, byte=isbyte
314 #endif
315 DO j = jMin, jMax
316 DO i = iMin, iMax
317 IF (frzmltMxL(i,j).GT.0. _d 0) THEN
318 C-- Net fluxes :
319 #ifdef ALLOW_ATM2D
320 pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
321 pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
322 sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
323 #else
324 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
325 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)
326 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
327 #endif
328
329 #ifdef ALLOW_DBUG_THSICE
330 IF (dBug(i,j,bi,bj)) WRITE(6,1010)
331 & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
332 & iceFrac(i,j), iceHeight(i,j,bi,bj),
333 & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
334 #endif
335 ENDIF
336
337 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
338 & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
339 & .AND. iceFrac(i,j) .LE.0. _d 0
340 IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
341 iceMask(i,j,bi,bj)=iceFrac(i,j)
342 IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
343 & snowAge(i,j,bi,bj) = 0. _d 0
344 ELSE
345 iceMask(i,j,bi,bj) = 0. _d 0
346 iceHeight(i,j,bi,bj)= 0. _d 0
347 snowHeight(i,j,bi,bj)=0. _d 0
348 snowAge(i,j,bi,bj) = 0. _d 0
349 Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
350 Tice1(i,j,bi,bj) = 0. _d 0
351 Tice2(i,j,bi,bj) = 0. _d 0
352 Qice1(i,j,bi,bj) = Lfresh
353 Qice2(i,j,bi,bj) = Lfresh
354 ENDIF
355 ENDDO
356 ENDDO
357
358 # ifdef ALLOW_AUTODIFF_TAMC
359 CADJ STORE snowHeight(:,:,bi,bj) =
360 CADJ & comlev1_bibj, key=iicekey, byte=isbyte
361 # endif
362 DO j = jMin, jMax
363 DO i = iMin, iMax
364 C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
365 sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
366 & + iceHeight(i,j,bi,bj)*rhoi
367 & )*iceMask(i,j,bi,bj)
368 ENDDO
369 ENDDO
370
371 IF ( thSIceAdvScheme.GT.0 ) THEN
372 C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
373 DO j = jMin, jMax
374 DO i = iMin, iMax
375 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
376 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
377 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)
378 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
379 ENDIF
380 ENDDO
381 ENDDO
382 ENDIF
383
384 #ifdef ALLOW_BULK_FORCE
385 IF ( useBulkForce ) THEN
386 CALL BULKF_FLUX_ADJUST(
387 I bi, bj, iMin, iMax, jMin, jMax,
388 I isIceFree, myTime, myIter, myThid )
389 ENDIF
390 #endif /* ALLOW_BULK_FORCE */
391
392 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
393 #endif /* ALLOW_THSICE */
394
395 RETURN
396 END

  ViewVC Help
Powered by ViewVC 1.1.22