/[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.20 - (hide annotations) (download)
Mon Apr 16 22:38:24 2007 UTC (17 years ago) by heimbach
Branch: MAIN
Changes since 1.19: +41 -1 lines
First set of modifs for TAF-ing thsice.

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

  ViewVC Help
Powered by ViewVC 1.1.22