/[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.21 - (show annotations) (download)
Tue Apr 17 23:42:33 2007 UTC (17 years ago) by heimbach
Branch: MAIN
Changes since 1.20: +18 -3 lines
2nd set of modifs for thsice adjoint.

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

  ViewVC Help
Powered by ViewVC 1.1.22