/[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.27 - (show annotations) (download)
Sat Oct 16 12:29:39 2010 UTC (13 years, 6 months ago) by heimbach
Branch: MAIN
Changes since 1.26: +10 -5 lines
Initial changes for adjoint of pkg/thsice

1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.26 2007/10/01 13:38:08 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 CADJ STORE iceheight(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
146 CADJ STORE icfrwatm(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
147 CADJ STORE qice1(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
148 CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
149 CADJ STORE snowheight(:,:,bi,bj) = comlev1_bibj,key=iicekey,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=iicekey,byte=isbyte
227 #endif
228
229 CALL THSICE_CALC_THICKN(
230 I bi, bj, siLo, siHi, sjLo, sjHi,
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 C-- Net fluxes :
244 DO j = jMin, jMax
245 DO i = iMin, iMax
246 #ifdef ALLOW_AUTODIFF_TAMC
247 ikey_1 = i
248 & + sNx*(j-1)
249 & + sNx*sNy*act1
250 & + sNx*sNy*max1*act2
251 & + sNx*sNy*max1*max2*act3
252 & + sNx*sNy*max1*max2*max3*act4
253 #endif /* ALLOW_AUTODIFF_TAMC */
254 C--
255 #ifdef ALLOW_AUTODIFF_TAMC
256 CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
257 #endif
258 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
259 C- weighted average net fluxes:
260 #ifdef ALLOW_AUTODIFF_TAMC
261 CADJ STORE fsalt(i,j) = comlev1_thsice_1, key=ikey_1
262 CADJ STORE flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
263 CADJ STORE frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
264 CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
265 #endif
266 icFrac = iceMask(i,j,bi,bj)
267 opFrac= 1. _d 0-icFrac
268 #ifdef ALLOW_ATM2D
269 pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
270 pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
271 sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
272 #else
273 icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
274 & - opFrac*Qnet(i,j,bi,bj)
275 icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
276 & + opFrac*EmPmR(i,j,bi,bj)
277 Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
278 EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)
279 & + opFrac*EmPmR(i,j,bi,bj)
280 saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
281 #endif
282
283 #ifdef ALLOW_DBUG_THSICE
284 IF (dBug(i,j,bi,bj)) WRITE(6,1010)
285 & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
286 & iceFrac(i,j), iceHeight(i,j,bi,bj),
287 & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
288 #endif
289
290 ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN
291 icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
292 icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj)
293 ELSE
294 icFlxAtm(i,j,bi,bj) = 0. _d 0
295 icFrwAtm(i,j,bi,bj) = 0. _d 0
296 ENDIF
297 ENDDO
298 ENDDO
299
300 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
301 C part.3 : freezing of sea-water
302 C over ice-free fraction and what is left from ice-covered fraction
303 C-------
304 CALL THSICE_EXTEND(
305 I bi, bj, siLo, siHi, sjLo, sjHi,
306 I iMin,iMax, jMin,jMax, dBugFlag,
307 I frzmltMxL, tFrzOce,
308 I tOceMxL(siLo,sjLo,bi,bj),
309 U iceFrac, iceHeight(siLo,sjLo,bi,bj),
310 U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
311 U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
312 U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
313 O flx2oc, frw2oc, fsalt,
314 I myTime, myIter, myThid )
315
316 #ifdef ALLOW_AUTODIFF_TAMC
317 CADJ STORE snowHeight(:,:,bi,bj) =
318 CADJ & comlev1_bibj,key=iicekey,byte=isbyte
319 #endif
320 DO j = jMin, jMax
321 DO i = iMin, iMax
322 IF (frzmltMxL(i,j).GT.0. _d 0) THEN
323 C-- Net fluxes :
324 #ifdef ALLOW_ATM2D
325 pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
326 pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
327 sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
328 #else
329 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
330 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)
331 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
332 #endif
333
334 #ifdef ALLOW_DBUG_THSICE
335 IF (dBug(i,j,bi,bj)) WRITE(6,1010)
336 & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
337 & iceFrac(i,j), iceHeight(i,j,bi,bj),
338 & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
339 #endif
340 ENDIF
341
342 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
343 & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
344 & .AND. iceFrac(i,j) .LE.0. _d 0
345 IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
346 iceMask(i,j,bi,bj)=iceFrac(i,j)
347 IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
348 & snowAge(i,j,bi,bj) = 0. _d 0
349 ELSE
350 iceMask(i,j,bi,bj) = 0. _d 0
351 iceHeight(i,j,bi,bj)= 0. _d 0
352 snowHeight(i,j,bi,bj)=0. _d 0
353 snowAge(i,j,bi,bj) = 0. _d 0
354 Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
355 Tice1(i,j,bi,bj) = 0. _d 0
356 Tice2(i,j,bi,bj) = 0. _d 0
357 Qice1(i,j,bi,bj) = Lfresh
358 Qice2(i,j,bi,bj) = Lfresh
359 ENDIF
360 ENDDO
361 ENDDO
362
363 # ifdef ALLOW_AUTODIFF_TAMC
364 CADJ STORE snowHeight(:,:,bi,bj) =
365 CADJ & comlev1_bibj,key=iicekey,byte=isbyte
366 # endif
367 DO j = jMin, jMax
368 DO i = iMin, iMax
369 C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
370 sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
371 & + iceHeight(i,j,bi,bj)*rhoi
372 & )*iceMask(i,j,bi,bj)
373 ENDDO
374 ENDDO
375
376 IF ( thSIceAdvScheme.GT.0 ) THEN
377 C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
378 DO j = jMin, jMax
379 DO i = iMin, iMax
380 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
381 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
382 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)
383 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
384 ENDIF
385 ENDDO
386 ENDDO
387 ENDIF
388
389 #ifdef ALLOW_BULK_FORCE
390 IF ( useBulkForce ) THEN
391 CALL BULKF_FLUX_ADJUST(
392 I bi, bj, iMin, iMax, jMin, jMax,
393 I isIceFree, myTime, myIter, myThid )
394 ENDIF
395 #endif /* ALLOW_BULK_FORCE */
396
397 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
398 #endif /* ALLOW_THSICE */
399
400 RETURN
401 END

  ViewVC Help
Powered by ViewVC 1.1.22