/[MITgcm]/MITgcm_contrib/jscott/code_changed/thsice_step_fwd.F
ViewVC logotype

Annotation of /MITgcm_contrib/jscott/code_changed/thsice_step_fwd.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Fri Aug 11 19:29:23 2006 UTC (20 years ago) by jscott
Branch: MAIN
atm2d package

1 jscott 1.1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.17 2006/05/25 18:03:25 jmc Exp $
2     C $Name: $
3    
4     #include "THSICE_OPTIONS.h"
5     #ifdef ALLOW_ATM2D
6     # include "ctrparam.h"
7     #endif
8    
9     CBOP
10     C !ROUTINE: THSICE_STEP_FWD
11     C !INTERFACE:
12     SUBROUTINE THSICE_STEP_FWD(
13     I bi, bj, iMin, iMax, jMin, jMax,
14     I prcAtm,
15     I myTime, myIter, myThid )
16     C !DESCRIPTION: \bv
17     C *==========================================================*
18     C | S/R THSICE_STEP_FWD
19     C | o Step Forward Therm-SeaIce model.
20     C *==========================================================*
21     C \ev
22    
23     C !USES:
24     IMPLICIT NONE
25    
26     C === Global variables ===
27     #include "SIZE.h"
28     #include "EEPARAMS.h"
29     #include "PARAMS.h"
30     #include "FFIELDS.h"
31     #ifdef ALLOW_ATM2D
32     # include "ATMSIZE.h"
33     # include "ATM2D_VARS.h"
34     #endif
35     #include "THSICE_SIZE.h"
36     #include "THSICE_PARAMS.h"
37     #include "THSICE_VARS.h"
38     #include "THSICE_TAVE.h"
39     INTEGER siLo, siHi, sjLo, sjHi
40     PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
41     PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
42    
43     C !INPUT/OUTPUT PARAMETERS:
44     C === Routine arguments ===
45     C- input:
46     C bi,bj :: tile indices
47     C iMin,iMax :: computation domain: 1rst index range
48     C jMin,jMax :: computation domain: 2nd index range
49     C prcAtm :: total precip from the atmosphere [kg/m2/s]
50     C myTime :: current Time of simulation [s]
51     C myIter :: current Iteration number in simulation
52     C myThid :: my Thread Id number
53     C-- Use fluxes hold in commom blocks
54     C- input:
55     C icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
56     C icFlxAtm :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
57     C icFrwAtm :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
58     C- output
59     C icFlxAtm :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
60     C icFrwAtm :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
61     INTEGER bi,bj
62     INTEGER iMin, iMax
63     INTEGER jMin, jMax
64     _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     _RL myTime
66     INTEGER myIter
67     INTEGER myThid
68     CEOP
69    
70     #ifdef ALLOW_THSICE
71     C !LOCAL VARIABLES:
72     C === Local variables ===
73     C iceFrac :: fraction of grid area covered in ice
74     C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2]
75     C frw2oc :: fresh-water flux from the ice to the ocean
76     C fsalt :: mass salt flux to the ocean
77     C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
78     C tFrzOce :: sea-water freezing temperature [oC] (function of S)
79     C isIceFree :: true for ice-free grid-cell that remains ice-free
80     C ageFac :: snow aging factor [1]
81     C snowFac :: snowing refreshing-age factor [units of 1/snowPr]
82     LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83     _RL iceFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84     _RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85     _RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86     _RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87     _RL tFrzOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88     _RL frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89     _RL ageFac
90     _RL snowFac
91     _RL cphm
92     _RL opFrac, icFrac
93     #ifdef ALLOW_DIAGNOSTICS
94     _RL tmpFac
95     #endif
96     INTEGER i,j
97     LOGICAL dBugFlag
98    
99     C- define grid-point location where to print debugging values
100     #include "THSICE_DEBUG.h"
101    
102     1010 FORMAT(A,1P4E14.6)
103    
104     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105    
106     C- Initialise
107     PRINT *,'*** Top of thsice_step_fwd'
108     dBugFlag = debugLevel.GE.debLevB
109     DO j = 1-OLy, sNy+OLy
110     DO i = 1-OLx, sNx+OLx
111     isIceFree(i,j) = .FALSE.
112     #ifdef ALLOW_ATM2D
113     sFluxFromIce(i,j) = 0. _d 0
114     #else
115     saltFlux(i,j,bi,bj) = 0. _d 0
116     #endif
117     #ifdef ALLOW_AUTODIFF_TAMC
118     iceFrac(i,j) = 0.
119     #endif
120     ENDDO
121     ENDDO
122    
123     ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
124     snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
125     DO j = jMin, jMax
126     DO i = iMin, iMax
127     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
128     C-- Snow aging :
129     snowAge(i,j,bi,bj) = thSIce_deltaT
130     & + snowAge(i,j,bi,bj)*ageFac
131     IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
132     & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
133     & * EXP( - snowFac*snowPrc(i,j,bi,bj) )
134     c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
135     c & /hNewSnowAge )
136     C-------
137     C note: Any flux of mass (here fresh water) that enter or leave the system
138     C with a non zero energy HAS TO be counted: add snow precip.
139     icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
140     & - Lfresh*snowPrc(i,j,bi,bj)
141     C--
142     ENDIF
143     ENDDO
144     ENDDO
145    
146     #ifdef ALLOW_DIAGNOSTICS
147     IF ( useDiagnostics ) THEN
148     tmpFac = 1. _d 0
149     CALL DIAGNOSTICS_FRACT_FILL(
150     I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
151     I 0,1,1,bi,bj,myThid)
152     CALL DIAGNOSTICS_FRACT_FILL(
153     I siceAlb, iceMask,tmpFac,1,'SIalbedo',
154     I 0,1,1,bi,bj,myThid)
155     ENDIF
156     #endif /* ALLOW_DIAGNOSTICS */
157     #ifndef ALLOW_ATM2D
158     DO j = jMin, jMax
159     DO i = iMin, iMax
160     siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
161     ENDDO
162     ENDDO
163     #endif
164    
165     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
166     C part.2 : ice-covered fraction ;
167     C change in ice/snow thickness and ice-fraction
168     C note: can only reduce the ice-fraction but not increase it.
169     C-------
170     DO j = jMin, jMax
171     DO i = iMin, iMax
172    
173     tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
174     cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
175     frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
176     & * cphm/ocean_deltaT
177     iceFrac(i,j) = iceMask(i,j,bi,bj)
178     flx2oc(i,j) = icFlxSW(i,j,bi,bj)
179     if (frzmltMxL(i,j).GT.0.D0) THEN
180     if ((i.eq.4).and.(j.eq.43)) THEN
181     PRINT *,'in frzmlt loop at:',i,j
182     print *, frzmltMxL(i,j),iceFrac(i,j),
183     & flx2oc(i,j),tOceMxL(i,j,bi,bj)
184     endif
185     endif
186     C-------
187     #ifdef ALLOW_DBUG_THSICE
188     IF ( dBug(i,j,bi,bj) ) THEN
189     IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
190     WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
191     WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
192     & iceFrac(i,j), iceHeight(i,j,bi,bj),
193     & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
194     WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
195     & tOceMxL(i,j,bi,bj), tFrzOce(i,j),
196     & frzmltMxL(i,j), Qnet(i,j,bi,bj)
197     ENDIF
198     IF (iceFrac(i,j).GT.0.)
199     & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
200     & iceFrac(i,j), icFlxAtm(i,j,bi,bj),
201     & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
202     ENDIF
203     #endif
204     ENDDO
205     ENDDO
206    
207     CALL THSICE_CALC_THICKN(
208     I bi, bj, siLo, siHi, sjLo, sjHi,
209     I iMin,iMax, jMin,jMax, dBugFlag,
210     I iceMask(siLo,sjLo,bi,bj), tFrzOce,
211     I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
212     I snowPrc(siLo,sjLo,bi,bj), prcAtm,
213     I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
214     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
215     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
216     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
217     U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
218     O frw2oc, fsalt,
219     I myTime, myIter, myThid )
220    
221     C-- Net fluxes :
222     DO j = jMin, jMax
223     DO i = iMin, iMax
224     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
225     C- weighted average net fluxes:
226     icFrac = iceMask(i,j,bi,bj)
227     opFrac= 1. _d 0-icFrac
228     #ifdef ALLOW_ATM2D
229     pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
230     pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
231     sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
232     #else
233     icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
234     & - opFrac*Qnet(i,j,bi,bj)
235     icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
236     & + opFrac*rhofw*EmPmR(i,j,bi,bj)
237     Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
238     EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
239     & + opFrac*EmPmR(i,j,bi,bj)
240     saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
241     #endif
242    
243     #ifdef ALLOW_DBUG_THSICE
244     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
245     & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
246     & iceFrac(i,j), iceHeight(i,j,bi,bj),
247     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
248     #endif
249    
250     #ifndef ALLOW_ATM2D
251     ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
252     icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
253     icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
254     ELSE
255     icFlxAtm(i,j,bi,bj) = 0. _d 0
256     icFrwAtm(i,j,bi,bj) = 0. _d 0
257     #endif
258     ENDIF
259     ENDDO
260     ENDDO
261    
262     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
263     C part.3 : freezing of sea-water
264     C over ice-free fraction and what is left from ice-covered fraction
265     C-------
266     CALL THSICE_EXTEND(
267     I bi, bj, siLo, siHi, sjLo, sjHi,
268     I iMin,iMax, jMin,jMax, dBugFlag,
269     I frzmltMxL, tFrzOce,
270     I tOceMxL(siLo,sjLo,bi,bj),
271     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
272     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
273     U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
274     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
275     O flx2oc, frw2oc, fsalt,
276     I myTime, myIter, myThid )
277    
278     DO j = jMin, jMax
279     DO i = iMin, iMax
280     IF (frzmltMxL(i,j).GT.0. _d 0) THEN
281     C-- Net fluxes :
282     #ifdef ALLOW_ATM2D
283     pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
284     pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
285     sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
286     #else
287     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
288     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
289     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
290     #endif
291    
292     #ifdef ALLOW_DBUG_THSICE
293     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
294     & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
295     & iceFrac(i,j), iceHeight(i,j,bi,bj),
296     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
297     #endif
298     ENDIF
299    
300     IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
301     & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
302     & .AND. iceFrac(i,j) .LE.0. _d 0
303     IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
304     iceMask(i,j,bi,bj)=iceFrac(i,j)
305     IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
306     & snowAge(i,j,bi,bj) = 0. _d 0
307     ELSE
308     iceMask(i,j,bi,bj) = 0. _d 0
309     iceHeight(i,j,bi,bj)= 0. _d 0
310     snowHeight(i,j,bi,bj)=0. _d 0
311     snowAge(i,j,bi,bj) = 0. _d 0
312     Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
313     Tice1(i,j,bi,bj) = 0. _d 0
314     Tice2(i,j,bi,bj) = 0. _d 0
315     Qice1(i,j,bi,bj) = 0. _d 0
316     Qice2(i,j,bi,bj) = 0. _d 0
317     ENDIF
318    
319     #ifdef ATMOSPHERIC_LOADING
320     C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
321     #ifndef ALLOW_ATM2D
322     sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
323     & + iceHeight(i,j,bi,bj)*rhoi
324     & )*iceMask(i,j,bi,bj)
325     #endif
326     #endif
327    
328     ENDDO
329     ENDDO
330    
331     #ifdef ALLOW_BULK_FORCE
332     IF ( useBulkForce ) THEN
333     CALL BULKF_FLUX_ADJUST(
334     I bi, bj, iMin, iMax, jMin, jMax,
335     I isIceFree, myTime, myIter, myThid )
336     ENDIF
337     #endif /* ALLOW_BULK_FORCE */
338    
339     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
340     #endif /* ALLOW_THSICE */
341    
342     CALL THSICE_AVE(1,1, myTime, myIter, myThid )
343    
344     RETURN
345     END

  ViewVC Help
Powered by ViewVC 1.1.22