/[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.3 - (hide annotations) (download)
Wed Apr 7 23:40:34 2004 UTC (20 years, 1 month ago) by jmc
Branch: MAIN
Changes since 1.2: +204 -150 lines
major changes in pkg/thsice: allows atmospheric model (AIM) to use thsice.
- split thsice_therm.F in 2 S/R: thsice_solve4temp.F & thsice_calc_thickn.F
- move most of the ocean & bulk_force interface in thsice_main.F
- add a "slab ocean" component to be used with atmospheric model

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

  ViewVC Help
Powered by ViewVC 1.1.22