/[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.18 - (hide annotations) (download)
Wed Apr 4 02:40:42 2007 UTC (17 years, 1 month ago) by jmc
Branch: MAIN
Changes since 1.17: +16 -2 lines
code to advect pkg/thSIce fields (testing is in progress).

1 jmc 1.18 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.17 2006/05/25 18:03:25 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 jmc 1.17 INTEGER siLo, siHi, sjLo, sjHi
34     PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
35     PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
36 jmc 1.13
37 jmc 1.1 C !INPUT/OUTPUT PARAMETERS:
38     C === Routine arguments ===
39 jmc 1.17 C- input:
40 jmc 1.3 C bi,bj :: tile indices
41     C iMin,iMax :: computation domain: 1rst index range
42     C jMin,jMax :: computation domain: 2nd index range
43 jmc 1.17 C prcAtm :: total precip from the atmosphere [kg/m2/s]
44     C myTime :: current Time of simulation [s]
45     C myIter :: current Iteration number in simulation
46     C myThid :: my Thread Id number
47     C-- Use fluxes hold in commom blocks
48 jmc 1.3 C- input:
49 jmc 1.17 C icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
50     C icFlxAtm :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
51     C icFrwAtm :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
52 jmc 1.3 C- output
53 jmc 1.17 C icFlxAtm :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
54     C icFrwAtm :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
55 jmc 1.1 INTEGER bi,bj
56     INTEGER iMin, iMax
57     INTEGER jMin, jMax
58 jmc 1.3 _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59 jmc 1.1 _RL myTime
60     INTEGER myIter
61     INTEGER myThid
62 jmc 1.3 CEOP
63 jmc 1.1
64     #ifdef ALLOW_THSICE
65     C !LOCAL VARIABLES:
66     C === Local variables ===
67 jmc 1.17 C iceFrac :: fraction of grid area covered in ice
68 jmc 1.3 C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2]
69     C frw2oc :: fresh-water flux from the ice to the ocean
70     C fsalt :: mass salt flux to the ocean
71     C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
72 jmc 1.17 C tFrzOce :: sea-water freezing temperature [oC] (function of S)
73 jmc 1.14 C isIceFree :: true for ice-free grid-cell that remains ice-free
74 jmc 1.17 C ageFac :: snow aging factor [1]
75     C snowFac :: snowing refreshing-age factor [units of 1/snowPr]
76     LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77     _RL iceFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78     _RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79     _RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80     _RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81     _RL tFrzOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82     _RL frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83     _RL ageFac
84     _RL snowFac
85     _RL cphm
86 jmc 1.3 _RL opFrac, icFrac
87 jmc 1.9 #ifdef ALLOW_DIAGNOSTICS
88 jmc 1.10 _RL tmpFac
89 jmc 1.9 #endif
90 jmc 1.17 INTEGER i,j
91     LOGICAL dBugFlag
92 jmc 1.1
93 jmc 1.17 C- define grid-point location where to print debugging values
94     #include "THSICE_DEBUG.h"
95 jmc 1.1
96 jmc 1.13 1010 FORMAT(A,1P4E14.6)
97 jmc 1.17
98     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
99    
100     C- Initialise
101     dBugFlag = debugLevel.GE.debLevB
102     DO j = 1-OLy, sNy+OLy
103     DO i = 1-OLx, sNx+OLx
104 jmc 1.14 isIceFree(i,j) = .FALSE.
105 jmc 1.17 saltFlux(i,j,bi,bj) = 0. _d 0
106 heimbach 1.16 #ifdef ALLOW_AUTODIFF_TAMC
107 jmc 1.17 iceFrac(i,j) = 0.
108 heimbach 1.16 #endif
109 jmc 1.4 ENDDO
110 jmc 1.17 ENDDO
111 jmc 1.1
112 jmc 1.17 ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
113     snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
114     DO j = jMin, jMax
115     DO i = iMin, iMax
116     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
117     C-- Snow aging :
118     snowAge(i,j,bi,bj) = thSIce_deltaT
119     & + snowAge(i,j,bi,bj)*ageFac
120     IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
121     & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
122     & * EXP( - snowFac*snowPrc(i,j,bi,bj) )
123     c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
124     c & /hNewSnowAge )
125 jmc 1.3 C-------
126 jmc 1.17 C note: Any flux of mass (here fresh water) that enter or leave the system
127     C with a non zero energy HAS TO be counted: add snow precip.
128     icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
129     & - Lfresh*snowPrc(i,j,bi,bj)
130     C--
131     ENDIF
132 jmc 1.3 ENDDO
133 jmc 1.17 ENDDO
134 jmc 1.3
135 jmc 1.9 #ifdef ALLOW_DIAGNOSTICS
136     IF ( useDiagnostics ) THEN
137 jmc 1.10 tmpFac = 1. _d 0
138     CALL DIAGNOSTICS_FRACT_FILL(
139 jmc 1.11 I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
140 jmc 1.10 I 0,1,1,bi,bj,myThid)
141 jmc 1.12 CALL DIAGNOSTICS_FRACT_FILL(
142     I siceAlb, iceMask,tmpFac,1,'SIalbedo',
143     I 0,1,1,bi,bj,myThid)
144 jmc 1.9 ENDIF
145     #endif /* ALLOW_DIAGNOSTICS */
146 jmc 1.12 DO j = jMin, jMax
147     DO i = iMin, iMax
148     siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
149     ENDDO
150     ENDDO
151 jmc 1.9
152 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
153 jmc 1.13 C part.2 : ice-covered fraction ;
154     C change in ice/snow thickness and ice-fraction
155 jmc 1.3 C note: can only reduce the ice-fraction but not increase it.
156     C-------
157 jmc 1.1 DO j = jMin, jMax
158     DO i = iMin, iMax
159 jmc 1.3
160 jmc 1.17 tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
161 jmc 1.3 cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
162 jmc 1.17 frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
163     & * cphm/ocean_deltaT
164     iceFrac(i,j) = iceMask(i,j,bi,bj)
165     flx2oc(i,j) = icFlxSW(i,j,bi,bj)
166 jmc 1.3 C-------
167 jmc 1.17 #ifdef ALLOW_DBUG_THSICE
168     IF ( dBug(i,j,bi,bj) ) THEN
169     IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
170 jmc 1.5 WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
171     WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
172 jmc 1.17 & iceFrac(i,j), iceHeight(i,j,bi,bj),
173 jmc 1.7 & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
174 jmc 1.17 WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
175     & tOceMxL(i,j,bi,bj), tFrzOce(i,j),
176     & frzmltMxL(i,j), Qnet(i,j,bi,bj)
177     ENDIF
178     IF (iceFrac(i,j).GT.0.)
179     & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
180     & iceFrac(i,j), icFlxAtm(i,j,bi,bj),
181     & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
182 jmc 1.1 ENDIF
183 jmc 1.17 #endif
184     ENDDO
185     ENDDO
186 jmc 1.1
187 jmc 1.17 CALL THSICE_CALC_THICKN(
188     I bi, bj, siLo, siHi, sjLo, sjHi,
189     I iMin,iMax, jMin,jMax, dBugFlag,
190     I iceMask(siLo,sjLo,bi,bj), tFrzOce,
191     I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
192     I snowPrc(siLo,sjLo,bi,bj), prcAtm,
193     I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
194     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
195     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
196     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
197     U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
198     O frw2oc, fsalt,
199     I myTime, myIter, myThid )
200 jmc 1.1
201     C-- Net fluxes :
202 jmc 1.17 DO j = jMin, jMax
203     DO i = iMin, iMax
204     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
205 jmc 1.3 C- weighted average net fluxes:
206     icFrac = iceMask(i,j,bi,bj)
207     opFrac= 1. _d 0-icFrac
208 jmc 1.17 icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
209     & - opFrac*Qnet(i,j,bi,bj)
210     icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
211     & + opFrac*rhofw*EmPmR(i,j,bi,bj)
212     Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
213     EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
214     & + opFrac*EmPmR(i,j,bi,bj)
215     saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
216    
217     #ifdef ALLOW_DBUG_THSICE
218     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
219     & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
220     & iceFrac(i,j), iceHeight(i,j,bi,bj),
221     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
222     #endif
223 jmc 1.1
224 jmc 1.3 ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
225 jmc 1.17 icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
226     icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
227 jmc 1.3 ELSE
228 jmc 1.17 icFlxAtm(i,j,bi,bj) = 0. _d 0
229     icFrwAtm(i,j,bi,bj) = 0. _d 0
230 jmc 1.1 ENDIF
231 jmc 1.17 ENDDO
232     ENDDO
233 jmc 1.1
234     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
235 jmc 1.3 C part.3 : freezing of sea-water
236 jmc 1.1 C over ice-free fraction and what is left from ice-covered fraction
237     C-------
238 jmc 1.17 CALL THSICE_EXTEND(
239     I bi, bj, siLo, siHi, sjLo, sjHi,
240     I iMin,iMax, jMin,jMax, dBugFlag,
241     I frzmltMxL, tFrzOce,
242     I tOceMxL(siLo,sjLo,bi,bj),
243     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
244     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
245     U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
246     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
247     O flx2oc, frw2oc, fsalt,
248     I myTime, myIter, myThid )
249    
250     DO j = jMin, jMax
251     DO i = iMin, iMax
252     IF (frzmltMxL(i,j).GT.0. _d 0) THEN
253 jmc 1.1 C-- Net fluxes :
254 jmc 1.17 Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
255     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
256     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
257    
258     #ifdef ALLOW_DBUG_THSICE
259     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
260     & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
261     & iceFrac(i,j), iceHeight(i,j,bi,bj),
262     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
263     #endif
264 jmc 1.1 ENDIF
265    
266 jmc 1.14 IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
267     & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
268 jmc 1.18 & .AND. iceFrac(i,j) .LE.0. _d 0
269 jmc 1.17 IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
270     iceMask(i,j,bi,bj)=iceFrac(i,j)
271     IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
272     & snowAge(i,j,bi,bj) = 0. _d 0
273 jmc 1.1 ELSE
274     iceMask(i,j,bi,bj) = 0. _d 0
275     iceHeight(i,j,bi,bj)= 0. _d 0
276     snowHeight(i,j,bi,bj)=0. _d 0
277 jmc 1.3 snowAge(i,j,bi,bj) = 0. _d 0
278 jmc 1.17 Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
279 jmc 1.1 Tice1(i,j,bi,bj) = 0. _d 0
280     Tice2(i,j,bi,bj) = 0. _d 0
281     Qice1(i,j,bi,bj) = 0. _d 0
282     Qice2(i,j,bi,bj) = 0. _d 0
283     ENDIF
284    
285 jmc 1.6 #ifdef ATMOSPHERIC_LOADING
286     C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
287     sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
288     & + iceHeight(i,j,bi,bj)*rhoi
289     & )*iceMask(i,j,bi,bj)
290     #endif
291    
292 jmc 1.1 ENDDO
293     ENDDO
294    
295 jmc 1.18 IF ( thSIceAdvScheme.GT.0 ) THEN
296     C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
297     DO j = jMin, jMax
298     DO i = iMin, iMax
299     IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
300     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
301     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)/rhofw
302     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
303     ENDIF
304     ENDDO
305     ENDDO
306     ENDIF
307    
308 jmc 1.14 #ifdef ALLOW_BULK_FORCE
309     IF ( useBulkForce ) THEN
310     CALL BULKF_FLUX_ADJUST(
311     I bi, bj, iMin, iMax, jMin, jMax,
312     I isIceFree, myTime, myIter, myThid )
313     ENDIF
314     #endif /* ALLOW_BULK_FORCE */
315    
316 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
317     #endif /* ALLOW_THSICE */
318    
319     RETURN
320     END

  ViewVC Help
Powered by ViewVC 1.1.22