/[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.5 - (hide annotations) (download)
Fri May 7 21:33:34 2004 UTC (20 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint53d_post, checkpoint54a_pre, checkpoint54a_post, checkpoint53c_post, checkpoint53b_post, checkpoint54b_post, checkpoint53b_pre, checkpoint53a_post, checkpoint54, checkpoint53g_post, checkpoint53f_post, checkpoint53d_pre
Changes since 1.4: +30 -14 lines
improve "debug" writing.

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

  ViewVC Help
Powered by ViewVC 1.1.22