/[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.21 - (hide 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 heimbach 1.21 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.20 2007/04/16 22:38:24 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     CALL DIAGNOSTICS_FRACT_FILL(
161 jmc 1.11 I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
162 jmc 1.10 I 0,1,1,bi,bj,myThid)
163 jmc 1.12 CALL DIAGNOSTICS_FRACT_FILL(
164     I siceAlb, iceMask,tmpFac,1,'SIalbedo',
165     I 0,1,1,bi,bj,myThid)
166 jmc 1.9 ENDIF
167     #endif /* ALLOW_DIAGNOSTICS */
168 jmc 1.12 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 jmc 1.9
174 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
175 jmc 1.13 C part.2 : ice-covered fraction ;
176     C change in ice/snow thickness and ice-fraction
177 jmc 1.3 C note: can only reduce the ice-fraction but not increase it.
178     C-------
179 jmc 1.1 DO j = jMin, jMax
180     DO i = iMin, iMax
181 jmc 1.3
182 jmc 1.17 tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
183 jmc 1.3 cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
184 jmc 1.17 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 jmc 1.3 C-------
189 jmc 1.17 #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 jmc 1.5 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 jmc 1.17 & iceFrac(i,j), iceHeight(i,j,bi,bj),
195 jmc 1.7 & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
196 jmc 1.17 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 jmc 1.1 ENDIF
205 jmc 1.17 #endif
206     ENDDO
207     ENDDO
208 jmc 1.1
209 heimbach 1.20 #ifdef ALLOW_AUTODIFF_TAMC
210     CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
211     #endif
212    
213 jmc 1.17 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 jmc 1.1
227     C-- Net fluxes :
228 jmc 1.17 DO j = jMin, jMax
229     DO i = iMin, iMax
230 heimbach 1.20 #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 jmc 1.17 IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
243 jmc 1.3 C- weighted average net fluxes:
244 heimbach 1.20 #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 jmc 1.3 icFrac = iceMask(i,j,bi,bj)
251     opFrac= 1. _d 0-icFrac
252 jmc 1.17 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 jmc 1.1
268 jmc 1.3 ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
269 jmc 1.17 icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
270     icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
271 jmc 1.3 ELSE
272 jmc 1.17 icFlxAtm(i,j,bi,bj) = 0. _d 0
273     icFrwAtm(i,j,bi,bj) = 0. _d 0
274 jmc 1.1 ENDIF
275 jmc 1.17 ENDDO
276     ENDDO
277 jmc 1.1
278     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
279 jmc 1.3 C part.3 : freezing of sea-water
280 jmc 1.1 C over ice-free fraction and what is left from ice-covered fraction
281     C-------
282 jmc 1.17 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 heimbach 1.21 #ifdef ALLOW_AUTODIFF_TAMC
295     CADJ STORE snowHeight(:,:,bi,bj) =
296     CADJ & comlev1_bibj, key=iicekey, byte=isbyte
297     #endif
298 jmc 1.17 DO j = jMin, jMax
299     DO i = iMin, iMax
300     IF (frzmltMxL(i,j).GT.0. _d 0) THEN
301 jmc 1.1 C-- Net fluxes :
302 jmc 1.17 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 jmc 1.1 ENDIF
313    
314 jmc 1.14 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
315     & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
316 jmc 1.18 & .AND. iceFrac(i,j) .LE.0. _d 0
317 jmc 1.17 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 jmc 1.1 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 jmc 1.3 snowAge(i,j,bi,bj) = 0. _d 0
326 jmc 1.17 Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
327 jmc 1.1 Tice1(i,j,bi,bj) = 0. _d 0
328     Tice2(i,j,bi,bj) = 0. _d 0
329 jmc 1.19 Qice1(i,j,bi,bj) = Lfresh
330     Qice2(i,j,bi,bj) = Lfresh
331 jmc 1.1 ENDIF
332 heimbach 1.21 ENDDO
333     ENDDO
334 jmc 1.1
335 jmc 1.6 #ifdef ATMOSPHERIC_LOADING
336 heimbach 1.21 # 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 jmc 1.6 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 jmc 1.1 ENDDO
347     ENDDO
348 heimbach 1.21 #endif
349 jmc 1.1
350 jmc 1.18 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 jmc 1.14 #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 jmc 1.1 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