/[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.24 - (show annotations) (download)
Fri Apr 27 15:50:41 2007 UTC (17 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59e, checkpoint59d, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b
Changes since 1.23: +1 -3 lines
make pLoad & sIceLoad always available (but still only used if
      ATMOSPHERIC_LOADING is #define).

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

  ViewVC Help
Powered by ViewVC 1.1.22