/[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.16 - (show annotations) (download)
Sun Apr 9 17:35:30 2006 UTC (18 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint58e_post, checkpoint58d_post
Changes since 1.15: +39 -28 lines
Starting thsice adjoint

1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.15 2006/03/26 00:45:40 jmc 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 U evpAtm, flxSW,
13 I myTime, myIter, myThid )
14 C !DESCRIPTION: \bv
15 C *==========================================================*
16 C | S/R THSICE_STEP_FWD
17 C | o Step Forward Therm-SeaIce model.
18 C *==========================================================*
19 C \ev
20
21 C !USES:
22 IMPLICIT NONE
23
24 C === Global variables ===
25 #include "SIZE.h"
26 #include "EEPARAMS.h"
27 #include "PARAMS.h"
28 #include "FFIELDS.h"
29 #include "THSICE_SIZE.h"
30 #include "THSICE_PARAMS.h"
31 #include "THSICE_VARS.h"
32 #include "THSICE_TAVE.h"
33
34 C !INPUT/OUTPUT PARAMETERS:
35 C === Routine arguments ===
36 C bi,bj :: tile indices
37 C iMin,iMax :: computation domain: 1rst index range
38 C jMin,jMax :: computation domain: 2nd index range
39 C- input:
40 C prcAtm :: total precip from the atmosphere [kg/m2/s]
41 C evpAtm :: (Inp) evaporation to the atmosphere [kg/m2/s] (>0 if evaporate)
42 C flxSW :: (Inp) short-wave heat flux (+=down): downward comp. only
43 C (part.1), becomes net SW flux into ocean (part.2).
44 C- output
45 C evpAtm :: (Out) net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
46 C flxSW :: (Out) net surf. heat flux from the atmosphere [W/m2], (+=down)
47 C myTime :: time counter for this thread
48 C myIter :: iteration counter for this thread
49 C myThid :: thread number for this instance of the routine.
50 INTEGER bi,bj
51 INTEGER iMin, iMax
52 INTEGER jMin, jMax
53 _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54 _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55 _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56 _RL myTime
57 INTEGER myIter
58 INTEGER myThid
59 CEOP
60
61 #ifdef ALLOW_THSICE
62 C !LOCAL VARIABLES:
63 C === Local variables ===
64 C snowPr :: snow precipitation [kg/m2/s]
65 C agingTime :: aging time scale (s)
66 C ageFac :: snow aging factor [1]
67 C albedo :: surface albedo [0-1]
68 C flxAtm :: net heat flux from the atmosphere (+=down) [W/m2]
69 C frwAtm :: net fresh-water flux (E-P) to the atmosphere [kg/m2/s]
70 C Fbot :: the oceanic heat flux already incorporated (ice_therm)
71 C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2]
72 C frw2oc :: fresh-water flux from the ice to the ocean
73 C fsalt :: mass salt flux to the ocean
74 C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
75 C TFrzOce :: sea-water freezing temperature [oC] (function of S)
76 C isIceFree :: true for ice-free grid-cell that remains ice-free
77 INTEGER i,j
78 _RL snowPr
79 _RL agingTime, ageFac
80 _RL albedo
81 _RL frwAtm
82 _RL flx2oc
83 _RL frw2oc
84 _RL fsalt
85 _RL TFrzOce, cphm, frzmltMxL
86 _RL Fbot, esurp
87 _RL opFrac, icFrac
88 _RL oceV2s, oceTs
89 _RL Tsf, Tice(nlyr), qicen(nlyr)
90 _RL tmpflx(0:2), tmpdTs
91 _RL hSnow
92 c
93 _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94 _RL compact(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95 _RL hIce(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96 c
97 LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98 #ifdef ALLOW_DIAGNOSTICS
99 _RL tmpFac
100 #endif
101
102 LOGICAL dBug
103
104 1010 FORMAT(A,1P4E14.6)
105 dBug = .FALSE.
106 C- Initialise flxAtm
107 DO j = 1-Oly, sNy+Oly
108 DO i = 1-Olx, sNx+Olx
109 flxAtm(i,j) = 0.
110 isIceFree(i,j) = .FALSE.
111 #ifdef ALLOW_AUTODIFF_TAMC
112 compact(i,j) = 0.
113 hIce(i,j) = 0.
114 #endif
115 ENDDO
116 ENDDO
117
118 IF ( fluidIsWater ) THEN
119 DO j = jMin, jMax
120 DO i = iMin, iMax
121 c dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )
122
123 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
124 C part.1 : ice-covered fraction ;
125 C Solve for surface and ice temperature (implicitly) ; compute surf. fluxes
126 C-------
127 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
128 icFrac = iceMask(i,j,bi,bj)
129 TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
130 hIce(i,j) = iceHeight(i,j,bi,bj)
131 hSnow = snowHeight(i,j,bi,bj)
132 Tsf = Tsrf(i,j,bi,bj)
133 qicen(1)= Qice1(i,j,bi,bj)
134 qicen(2)= Qice2(i,j,bi,bj)
135 IF ( dBug ) THEN
136 WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
137 WRITE(6,1010) 'ThSI_FWD:-0- iceMask, hIc, hSn, Tsf =',
138 & icFrac, hIce,hSnow,Tsf
139 ENDIF
140
141 CALL THSICE_ALBEDO(
142 I hIce(i,j), hSnow, Tsf, snowAge(i,j,bi,bj),
143 O albedo,
144 I myThid )
145 flxSW(i,j) = flxSW(i,j)*(1. _d 0 - albedo)
146 siceAlb(i,j,bi,bj) = albedo
147
148 CALL THSICE_SOLVE4TEMP(
149 I useBulkForce, tmpflx, TFrzOce, hIce(i,j), hSnow,
150 U flxSW(i,j), Tsf, qicen,
151 O Tice, sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
152 O tmpdTs, flxAtm(i,j), evpAtm(i,j),
153 I i,j, bi,bj, myThid)
154
155 #ifdef SHORTWAVE_HEATING
156 C-- Update Fluxes :
157 opFrac= 1. _d 0-icFrac
158 Qsw(i,j,bi,bj)=-icFrac*flxSW(i,j) +opFrac*Qsw(i,j,bi,bj)
159 #endif
160 C-- Update Sea-Ice state :
161 Tsrf(i,j,bi,bj) =Tsf
162 Tice1(i,j,bi,bj)=Tice(1)
163 Tice2(i,j,bi,bj)=Tice(2)
164 Qice1(i,j,bi,bj)=qicen(1)
165 Qice2(i,j,bi,bj)=qicen(2)
166 IF ( dBug ) THEN
167 WRITE(6,1010) 'ThSI_FWD: Tsf, Tice(1,2), TFrzOce =',
168 & Tsf, Tice, TFrzOce
169 WRITE(6,1010) 'ThSI_FWD: sHeat,fxCndBt, fxAtm,evAtm=',
170 & sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj),
171 & flxAtm(i,j), evpAtm(i,j)
172 ENDIF
173 ENDIF
174 ENDDO
175 ENDDO
176 ENDIF
177 dBug = .FALSE.
178
179 #ifdef ALLOW_DIAGNOSTICS
180 IF ( useDiagnostics ) THEN
181 tmpFac = 1. _d 0
182 CALL DIAGNOSTICS_FRACT_FILL(
183 I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
184 I 0,1,1,bi,bj,myThid)
185 CALL DIAGNOSTICS_FRACT_FILL(
186 I siceAlb, iceMask,tmpFac,1,'SIalbedo',
187 I 0,1,1,bi,bj,myThid)
188 ENDIF
189 #endif /* ALLOW_DIAGNOSTICS */
190 DO j = jMin, jMax
191 DO i = iMin, iMax
192 siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
193 ENDDO
194 ENDDO
195
196 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
197 C part.2 : ice-covered fraction ;
198 C change in ice/snow thickness and ice-fraction
199 C note: can only reduce the ice-fraction but not increase it.
200 C-------
201 agingTime = 50. _d 0 * 86400. _d 0
202 ageFac = 1. _d 0 - thSIce_deltaT/agingTime
203 DO j = jMin, jMax
204 DO i = iMin, iMax
205 c dBug = ( bi.EQ.3 .AND. i.EQ.15 .AND. j.EQ.11 )
206
207 TFrzOce = -mu_Tf*sOceMxL(i,j,bi,bj)
208 oceTs = tOceMxL(i,j,bi,bj)
209 cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
210 frzmltMxL = (TFrzOce-oceTs)*cphm/ocean_deltaT
211
212 Fbot = 0. _d 0
213 saltFlux(i,j,bi,bj) = 0. _d 0
214 compact(i,j)= iceMask(i,j,bi,bj)
215 C-------
216 IF (dBug .AND. (frzmltMxL.GT.0. .OR. compact(i,j).GT.0.) ) THEN
217 WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
218 WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
219 & compact(i,j), iceHeight(i,j,bi,bj),
220 & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
221 WRITE(6,1010) 'ThSI_FWD: ocTs,TFrzOce,frzmltMxL,Qnet=',
222 & oceTs, TFrzOce, frzmltMxL,Qnet(i,j,bi,bj)
223 ENDIF
224 C-------
225 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
226
227 oceV2s = v2ocMxL(i,j,bi,bj)
228 snowPr = snowPrc(i,j,bi,bj)
229 hIce(i,j) = iceHeight(i,j,bi,bj)
230 hSnow = snowHeight(i,j,bi,bj)
231 Tsf = Tsrf(i,j,bi,bj)
232 qicen(1)= Qice1(i,j,bi,bj)
233 qicen(2)= Qice2(i,j,bi,bj)
234 flx2oc = flxSW(i,j)
235
236 CALL THSICE_CALC_THICKN(
237 I frzmltMxL, TFrzOce, oceTs, oceV2s, snowPr,
238 I sHeating(i,j,bi,bj), flxCndBt(i,j,bi,bj), evpAtm(i,j),
239 U compact(i,j), hIce(i,j), hSnow, Tsf, qicen, flx2oc,
240 O frw2oc, fsalt, Fbot,
241 I dBug, myThid)
242
243 C- note : snowPr was not supposed to be modified in THSICE_THERM ;
244 C but to reproduce old results, is reset to zero if Tsf >= 0
245 snowPrc(i,j,bi,bj) = snowPr
246
247 C-- Snow aging :
248 snowAge(i,j,bi,bj) = thSIce_deltaT
249 & + snowAge(i,j,bi,bj)*ageFac
250 IF ( snowPr.GT.0. _d 0 )
251 & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
252 & * EXP( -(thSIce_deltaT*snowPr/rhos)/hNewSnowAge )
253 C--
254
255 C-- Diagnostic of Atmospheric Fluxes over sea-ice :
256 frwAtm = evpAtm(i,j) - prcAtm(i,j)
257 C note: Any flux of mass (here fresh water) that enter or leave the system
258 C with a non zero energy HAS TO be counted: add snow precip.
259 flxAtm(i,j) = flxAtm(i,j) - Lfresh*snowPrc(i,j,bi,bj)
260
261 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
262 IF (dBug) WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
263 & iceMask(i,j,bi,bj),flxAtm(i,j),evpAtm(i,j),-Lfresh*snowPr
264 IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc=',
265 & compact(i,j),flx2oc,fsalt,frw2oc
266 #ifdef CHECK_ENERGY_CONSERV
267 icFrac = iceMask(i,j,bi,bj)
268 CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 0,
269 I icFrac, compact(i,j), hIce(i,j), hSnow, qicen,
270 I flx2oc, frw2oc, fsalt, flxAtm(i,j), frwAtm,
271 I myTime, myIter, myThid )
272 #endif /* CHECK_ENERGY_CONSERV */
273 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
274
275 C-- Update Sea-Ice state :
276 c iceMask(i,j,bi,bj)=compact
277 iceHeight(i,j,bi,bj) = hIce(i,j)
278 snowHeight(i,j,bi,bj)= hSnow
279 Tsrf(i,j,bi,bj) =Tsf
280 Qice1(i,j,bi,bj)=qicen(1)
281 Qice2(i,j,bi,bj)=qicen(2)
282
283 C-- Net fluxes :
284 frw2oc = frw2oc + (prcAtm(i,j)-snowPrc(i,j,bi,bj))
285 C- weighted average net fluxes:
286 icFrac = iceMask(i,j,bi,bj)
287 opFrac= 1. _d 0-icFrac
288 flxAtm(i,j) = icFrac*flxAtm(i,j) - opFrac*Qnet(i,j,bi,bj)
289 frwAtm = icFrac*frwAtm + opFrac*rhofw*EmPmR(i,j,bi,bj)
290 Qnet(i,j,bi,bj)=-icFrac*flx2oc +opFrac*Qnet(i,j,bi,bj)
291 EmPmR(i,j,bi,bj)=-icFrac*frw2oc/rhofw+opFrac*EmPmR(i,j,bi,bj)
292 saltFlux(i,j,bi,bj)=-icFrac*fsalt
293
294 IF (dBug) WRITE(6,1010)
295 & 'ThSI_FWD:-3- compact, hIc, hSn, Qnet =',
296 & compact(i,j),hIce(i,j),hSnow,
297 & Qnet(i,j,bi,bj)
298
299 ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
300 flxAtm(i,j) = -Qnet(i,j,bi,bj)
301 frwAtm = rhofw*EmPmR(i,j,bi,bj)
302 ELSE
303 flxAtm(i,j) = 0. _d 0
304 frwAtm = 0. _d 0
305 ENDIF
306
307 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
308 C part.3 : freezing of sea-water
309 C over ice-free fraction and what is left from ice-covered fraction
310 C-------
311 c compact(i,j) = iceMask(i,j,bi,bj)
312 hIce(i,j) = iceHeight(i,j,bi,bj)
313 hSnow = snowHeight(i,j,bi,bj)
314
315 esurp = frzmltMxL - Fbot*iceMask(i,j,bi,bj)
316 IF (esurp.GT.0. _d 0) THEN
317 icFrac = compact(i,j)
318 qicen(1)= Qice1(i,j,bi,bj)
319 qicen(2)= Qice2(i,j,bi,bj)
320 CALL THSICE_EXTEND(
321 I esurp, TFrzOce,
322 U oceTs, compact(i,j), hIce(i,j), hSnow, qicen,
323 O flx2oc, frw2oc, fsalt,
324 I dBug, myThid )
325 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
326 IF (dBug) WRITE(6,1010) 'ThSI_FWD: compact,flx2oc,fsalt,frw2oc='
327 & ,compact(i,j),flx2oc,fsalt,frw2oc
328 #ifdef CHECK_ENERGY_CONSERV
329 tmpflx(1) = 0.
330 tmpflx(2) = 0.
331 CALL THSICE_CHECK_CONSERV( dBug, i, j, bi, bj, 1,
332 I icFrac, compact(i,j), hIce(i,j), hSnow, qicen,
333 I flx2oc, frw2oc, fsalt, tmpflx(1), tmpflx(2),
334 I myTime, myIter, myThid )
335 #endif /* CHECK_ENERGY_CONSERV */
336 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
337 C-- Update Sea-Ice state :
338 IF ( compact(i,j).GT.0. _d 0 .AND. icFrac.EQ.0. _d 0) THEN
339 Tsrf(i,j,bi,bj) = TFrzOce
340 Tice1(i,j,bi,bj) = TFrzOce
341 Tice2(i,j,bi,bj) = TFrzOce
342 Qice1(i,j,bi,bj) = qicen(1)
343 Qice2(i,j,bi,bj) = qicen(2)
344 ENDIF
345 iceHeight(i,j,bi,bj) = hIce(i,j)
346 snowHeight(i,j,bi,bj)= hSnow
347 C-- Net fluxes :
348 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc
349 EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc/rhofw
350 saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt
351
352 IF (dBug) WRITE(6,1010)
353 & 'ThSI_FWD:-4- compact, hIc, hSn, Qnet =',
354 & compact(i,j),hIce(i,j),hSnow,
355 & Qnet(i,j,bi,bj)
356 C-- - if esurp > 0 : end
357 ENDIF
358
359 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
360 & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
361 & .AND. compact(i,j) .LE.0. _d 0
362 IF ( compact(i,j) .GT. 0. _d 0 ) THEN
363 iceMask(i,j,bi,bj)=compact(i,j)
364 IF ( hSnow .EQ. 0. _d 0 ) snowAge(i,j,bi,bj) = 0. _d 0
365 ELSE
366 iceMask(i,j,bi,bj) = 0. _d 0
367 iceHeight(i,j,bi,bj)= 0. _d 0
368 snowHeight(i,j,bi,bj)=0. _d 0
369 snowAge(i,j,bi,bj) = 0. _d 0
370 Tsrf(i,j,bi,bj) = oceTs
371 Tice1(i,j,bi,bj) = 0. _d 0
372 Tice2(i,j,bi,bj) = 0. _d 0
373 Qice1(i,j,bi,bj) = 0. _d 0
374 Qice2(i,j,bi,bj) = 0. _d 0
375 ENDIF
376
377 C-- Return atmospheric fluxes in evpAtm & flxSW (same sign and units):
378 evpAtm(i,j) = frwAtm
379 flxSW (i,j) = flxAtm(i,j)
380
381 #ifdef ATMOSPHERIC_LOADING
382 C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
383 sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
384 & + iceHeight(i,j,bi,bj)*rhoi
385 & )*iceMask(i,j,bi,bj)
386 #endif
387
388 ENDDO
389 ENDDO
390
391 #ifdef ALLOW_BULK_FORCE
392 IF ( useBulkForce ) THEN
393 CALL BULKF_FLUX_ADJUST(
394 I bi, bj, iMin, iMax, jMin, jMax,
395 I isIceFree, myTime, myIter, myThid )
396 ENDIF
397 #endif /* ALLOW_BULK_FORCE */
398
399 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
400 #endif /* ALLOW_THSICE */
401
402 RETURN
403 END

  ViewVC Help
Powered by ViewVC 1.1.22