/[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.10 - (hide annotations) (download)
Sun Jul 10 01:24:25 2005 UTC (18 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.9: +6 -14 lines
use DIAGNOSTICS_FRACT_FILL insted of DIAGNOSTICS_FILL to get correct
 statistics diagnostics over sea-ice.

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

  ViewVC Help
Powered by ViewVC 1.1.22