/[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.15 - (hide annotations) (download)
Sun Mar 26 00:45:40 2006 UTC (18 years, 1 month ago) by jmc
Branch: MAIN
Changes since 1.14: +3 -3 lines
fix a debug write.

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

  ViewVC Help
Powered by ViewVC 1.1.22