/[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.9 - (hide annotations) (download)
Mon Jan 31 19:37:06 2005 UTC (19 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint57k_post, checkpoint57d_post, checkpoint57g_post, checkpoint57i_post, checkpoint57e_post, checkpoint57g_pre, checkpoint57f_pre, eckpoint57e_pre, checkpoint57h_done, checkpoint57f_post, checkpoint57j_post, checkpoint57h_pre, checkpoint57h_post
Changes since 1.8: +22 -5 lines
store sea-ice albedo in a common block (for diagnostics).

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

  ViewVC Help
Powered by ViewVC 1.1.22