/[MITgcm]/MITgcm/pkg/thsice/thsice_step_fwd.F
ViewVC logotype

Annotation of /MITgcm/pkg/thsice/thsice_step_fwd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.22 - (hide annotations) (download)
Thu Apr 19 04:31:28 2007 UTC (17 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59
Changes since 1.21: +2 -1 lines
fix albedo & snow-precip diagnostics (when using advection).

1 jmc 1.22 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.21 2007/04/17 23:42:33 heimbach Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "THSICE_OPTIONS.h"
5 jmc 1.3
6     CBOP
7 jmc 1.1 C !ROUTINE: THSICE_STEP_FWD
8     C !INTERFACE:
9 jmc 1.13 SUBROUTINE THSICE_STEP_FWD(
10 jmc 1.1 I bi, bj, iMin, iMax, jMin, jMax,
11 jmc 1.3 I prcAtm,
12 jmc 1.1 I myTime, myIter, myThid )
13 jmc 1.3 C !DESCRIPTION: \bv
14 jmc 1.1 C *==========================================================*
15 jmc 1.13 C | S/R THSICE_STEP_FWD
16 jmc 1.1 C | o Step Forward Therm-SeaIce model.
17     C *==========================================================*
18 jmc 1.3 C \ev
19 jmc 1.1
20     C !USES:
21     IMPLICIT NONE
22 jmc 1.3
23 jmc 1.1 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 jmc 1.3 #include "THSICE_VARS.h"
31     #include "THSICE_TAVE.h"
32 jmc 1.18 #include "THSICE_2DYN.h"
33 heimbach 1.20 #ifdef ALLOW_AUTODIFF_TAMC
34     # include "tamc.h"
35     # include "tamc_keys.h"
36     #endif
37    
38 jmc 1.17 INTEGER siLo, siHi, sjLo, sjHi
39     PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
40     PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
41 jmc 1.13
42 jmc 1.1 C !INPUT/OUTPUT PARAMETERS:
43     C === Routine arguments ===
44 jmc 1.17 C- input:
45 jmc 1.3 C bi,bj :: tile indices
46     C iMin,iMax :: computation domain: 1rst index range
47     C jMin,jMax :: computation domain: 2nd index range
48 jmc 1.17 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 jmc 1.3 C- input:
54 jmc 1.17 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 jmc 1.3 C- output
58 jmc 1.17 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 jmc 1.1 INTEGER bi,bj
61     INTEGER iMin, iMax
62     INTEGER jMin, jMax
63 jmc 1.3 _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 jmc 1.1 _RL myTime
65     INTEGER myIter
66     INTEGER myThid
67 jmc 1.3 CEOP
68 jmc 1.1
69     #ifdef ALLOW_THSICE
70     C !LOCAL VARIABLES:
71     C === Local variables ===
72 jmc 1.17 C iceFrac :: fraction of grid area covered in ice
73 jmc 1.3 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 jmc 1.17 C tFrzOce :: sea-water freezing temperature [oC] (function of S)
78 jmc 1.14 C isIceFree :: true for ice-free grid-cell that remains ice-free
79 jmc 1.17 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 jmc 1.3 _RL opFrac, icFrac
92 jmc 1.9 #ifdef ALLOW_DIAGNOSTICS
93 jmc 1.10 _RL tmpFac
94 jmc 1.9 #endif
95 jmc 1.17 INTEGER i,j
96     LOGICAL dBugFlag
97 jmc 1.1
98 jmc 1.17 C- define grid-point location where to print debugging values
99     #include "THSICE_DEBUG.h"
100 jmc 1.1
101 jmc 1.13 1010 FORMAT(A,1P4E14.6)
102 jmc 1.17
103     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
104    
105 heimbach 1.20 #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 jmc 1.17 C- Initialise
119     dBugFlag = debugLevel.GE.debLevB
120     DO j = 1-OLy, sNy+OLy
121     DO i = 1-OLx, sNx+OLx
122 jmc 1.14 isIceFree(i,j) = .FALSE.
123 jmc 1.17 saltFlux(i,j,bi,bj) = 0. _d 0
124 heimbach 1.16 #ifdef ALLOW_AUTODIFF_TAMC
125 jmc 1.17 iceFrac(i,j) = 0.
126 heimbach 1.16 #endif
127 jmc 1.4 ENDDO
128 jmc 1.17 ENDDO
129 jmc 1.1
130 jmc 1.17 ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
131     snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
132 heimbach 1.21
133     #ifdef ALLOW_AUTODIFF_TAMC
134     CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
135     #endif
136 jmc 1.17 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 jmc 1.3 C-------
148 jmc 1.17 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 jmc 1.3 ENDDO
155 jmc 1.17 ENDDO
156 jmc 1.3
157 jmc 1.9 #ifdef ALLOW_DIAGNOSTICS
158     IF ( useDiagnostics ) THEN
159 jmc 1.10 tmpFac = 1. _d 0
160 jmc 1.22 CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
161 jmc 1.10 CALL DIAGNOSTICS_FRACT_FILL(
162 jmc 1.11 I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
163 jmc 1.10 I 0,1,1,bi,bj,myThid)
164 jmc 1.12 CALL DIAGNOSTICS_FRACT_FILL(
165     I siceAlb, iceMask,tmpFac,1,'SIalbedo',
166     I 0,1,1,bi,bj,myThid)
167 jmc 1.9 ENDIF
168     #endif /* ALLOW_DIAGNOSTICS */
169 jmc 1.12 DO j = jMin, jMax
170     DO i = iMin, iMax
171     siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
172     ENDDO
173     ENDDO
174 jmc 1.9
175 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
176 jmc 1.13 C part.2 : ice-covered fraction ;
177     C change in ice/snow thickness and ice-fraction
178 jmc 1.3 C note: can only reduce the ice-fraction but not increase it.
179     C-------
180 jmc 1.1 DO j = jMin, jMax
181     DO i = iMin, iMax
182 jmc 1.3
183 jmc 1.17 tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
184 jmc 1.3 cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
185 jmc 1.17 frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
186     & * cphm/ocean_deltaT
187     iceFrac(i,j) = iceMask(i,j,bi,bj)
188     flx2oc(i,j) = icFlxSW(i,j,bi,bj)
189 jmc 1.3 C-------
190 jmc 1.17 #ifdef ALLOW_DBUG_THSICE
191     IF ( dBug(i,j,bi,bj) ) THEN
192     IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
193 jmc 1.5 WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
194     WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
195 jmc 1.17 & iceFrac(i,j), iceHeight(i,j,bi,bj),
196 jmc 1.7 & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
197 jmc 1.17 WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
198     & tOceMxL(i,j,bi,bj), tFrzOce(i,j),
199     & frzmltMxL(i,j), Qnet(i,j,bi,bj)
200     ENDIF
201     IF (iceFrac(i,j).GT.0.)
202     & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
203     & iceFrac(i,j), icFlxAtm(i,j,bi,bj),
204     & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
205 jmc 1.1 ENDIF
206 jmc 1.17 #endif
207     ENDDO
208     ENDDO
209 jmc 1.1
210 heimbach 1.20 #ifdef ALLOW_AUTODIFF_TAMC
211     CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
212     #endif
213    
214 jmc 1.17 CALL THSICE_CALC_THICKN(
215     I bi, bj, siLo, siHi, sjLo, sjHi,
216     I iMin,iMax, jMin,jMax, dBugFlag,
217     I iceMask(siLo,sjLo,bi,bj), tFrzOce,
218     I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
219     I snowPrc(siLo,sjLo,bi,bj), prcAtm,
220     I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
221     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
222     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
223     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
224     U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
225     O frw2oc, fsalt,
226     I myTime, myIter, myThid )
227 jmc 1.1
228     C-- Net fluxes :
229 jmc 1.17 DO j = jMin, jMax
230     DO i = iMin, iMax
231 heimbach 1.20 #ifdef ALLOW_AUTODIFF_TAMC
232     ikey_1 = i
233     & + sNx*(j-1)
234     & + sNx*sNy*act1
235     & + sNx*sNy*max1*act2
236     & + sNx*sNy*max1*max2*act3
237     & + sNx*sNy*max1*max2*max3*act4
238     #endif /* ALLOW_AUTODIFF_TAMC */
239     C--
240     #ifdef ALLOW_AUTODIFF_TAMC
241     CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
242     #endif
243 jmc 1.17 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
244 jmc 1.3 C- weighted average net fluxes:
245 heimbach 1.20 #ifdef ALLOW_AUTODIFF_TAMC
246     CADJ STORE fsalt(i,j) = comlev1_thsice_1, key=ikey_1
247     CADJ STORE flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
248     CADJ STORE frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
249     CADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
250     #endif
251 jmc 1.3 icFrac = iceMask(i,j,bi,bj)
252     opFrac= 1. _d 0-icFrac
253 jmc 1.17 icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
254     & - opFrac*Qnet(i,j,bi,bj)
255     icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
256     & + opFrac*rhofw*EmPmR(i,j,bi,bj)
257     Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
258     EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
259     & + opFrac*EmPmR(i,j,bi,bj)
260     saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
261    
262     #ifdef ALLOW_DBUG_THSICE
263     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
264     & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
265     & iceFrac(i,j), iceHeight(i,j,bi,bj),
266     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
267     #endif
268 jmc 1.1
269 jmc 1.3 ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
270 jmc 1.17 icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
271     icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
272 jmc 1.3 ELSE
273 jmc 1.17 icFlxAtm(i,j,bi,bj) = 0. _d 0
274     icFrwAtm(i,j,bi,bj) = 0. _d 0
275 jmc 1.1 ENDIF
276 jmc 1.17 ENDDO
277     ENDDO
278 jmc 1.1
279     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
280 jmc 1.3 C part.3 : freezing of sea-water
281 jmc 1.1 C over ice-free fraction and what is left from ice-covered fraction
282     C-------
283 jmc 1.17 CALL THSICE_EXTEND(
284     I bi, bj, siLo, siHi, sjLo, sjHi,
285     I iMin,iMax, jMin,jMax, dBugFlag,
286     I frzmltMxL, tFrzOce,
287     I tOceMxL(siLo,sjLo,bi,bj),
288     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
289     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
290     U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
291     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
292     O flx2oc, frw2oc, fsalt,
293     I myTime, myIter, myThid )
294    
295 heimbach 1.21 #ifdef ALLOW_AUTODIFF_TAMC
296     CADJ STORE snowHeight(:,:,bi,bj) =
297     CADJ & comlev1_bibj, key=iicekey, byte=isbyte
298     #endif
299 jmc 1.17 DO j = jMin, jMax
300     DO i = iMin, iMax
301     IF (frzmltMxL(i,j).GT.0. _d 0) THEN
302 jmc 1.1 C-- Net fluxes :
303 jmc 1.17 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
304     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
305     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
306    
307     #ifdef ALLOW_DBUG_THSICE
308     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
309     & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
310     & iceFrac(i,j), iceHeight(i,j,bi,bj),
311     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
312     #endif
313 jmc 1.1 ENDIF
314    
315 jmc 1.14 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
316     & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
317 jmc 1.18 & .AND. iceFrac(i,j) .LE.0. _d 0
318 jmc 1.17 IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
319     iceMask(i,j,bi,bj)=iceFrac(i,j)
320     IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
321     & snowAge(i,j,bi,bj) = 0. _d 0
322 jmc 1.1 ELSE
323     iceMask(i,j,bi,bj) = 0. _d 0
324     iceHeight(i,j,bi,bj)= 0. _d 0
325     snowHeight(i,j,bi,bj)=0. _d 0
326 jmc 1.3 snowAge(i,j,bi,bj) = 0. _d 0
327 jmc 1.17 Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
328 jmc 1.1 Tice1(i,j,bi,bj) = 0. _d 0
329     Tice2(i,j,bi,bj) = 0. _d 0
330 jmc 1.19 Qice1(i,j,bi,bj) = Lfresh
331     Qice2(i,j,bi,bj) = Lfresh
332 jmc 1.1 ENDIF
333 heimbach 1.21 ENDDO
334     ENDDO
335 jmc 1.1
336 jmc 1.6 #ifdef ATMOSPHERIC_LOADING
337 heimbach 1.21 # ifdef ALLOW_AUTODIFF_TAMC
338     CADJ STORE snowHeight(:,:,bi,bj) =
339     CADJ & comlev1_bibj, key=iicekey, byte=isbyte
340     # endif
341     DO j = jMin, jMax
342     DO i = iMin, iMax
343 jmc 1.6 C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
344     sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
345     & + iceHeight(i,j,bi,bj)*rhoi
346     & )*iceMask(i,j,bi,bj)
347 jmc 1.1 ENDDO
348     ENDDO
349 heimbach 1.21 #endif
350 jmc 1.1
351 jmc 1.18 IF ( thSIceAdvScheme.GT.0 ) THEN
352     C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
353     DO j = jMin, jMax
354     DO i = iMin, iMax
355     IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
356     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
357     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw
358     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
359     ENDIF
360     ENDDO
361     ENDDO
362     ENDIF
363    
364 jmc 1.14 #ifdef ALLOW_BULK_FORCE
365     IF ( useBulkForce ) THEN
366     CALL BULKF_FLUX_ADJUST(
367     I bi, bj, iMin, iMax, jMin, jMax,
368     I isIceFree, myTime, myIter, myThid )
369     ENDIF
370     #endif /* ALLOW_BULK_FORCE */
371    
372 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
373     #endif /* ALLOW_THSICE */
374    
375     RETURN
376     END

  ViewVC Help
Powered by ViewVC 1.1.22