/[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.4 - (hide annotations) (download)
Thu Apr 8 18:54:26 2004 UTC (20 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint52n_post, checkpoint52m_post, checkpoint53
Changes since 1.3: +8 -2 lines
disable inaccurate diagnostics.

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

  ViewVC Help
Powered by ViewVC 1.1.22