/[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.35 - (show annotations) (download)
Wed Feb 8 22:32:24 2012 UTC (12 years, 4 months ago) by jmc
Branch: MAIN
Changes since 1.34: +16 -20 lines
- diagnose directly amount of seawater freezing in THSICE_CALC_THICKN (not used).
- fix and simplify interface to salt_plume pkg.

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

  ViewVC Help
Powered by ViewVC 1.1.22