/[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.2 - (hide annotations) (download)
Tue Aug 22 20:28:31 2006 UTC (20 years ago) by jscott
Branch: MAIN
Changes since 1.1: +1 -4 lines
removed extra print statements

1 jscott 1.2 C $Header: /u/gcmpack/MITgcm_contrib/jscott/code_changed/thsice_step_fwd.F,v 1.1 2006/08/11 19:29:23 jscott Exp $
2 jscott 1.1 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     dBugFlag = debugLevel.GE.debLevB
108     DO j = 1-OLy, sNy+OLy
109     DO i = 1-OLx, sNx+OLx
110     isIceFree(i,j) = .FALSE.
111     #ifdef ALLOW_ATM2D
112     sFluxFromIce(i,j) = 0. _d 0
113     #else
114     saltFlux(i,j,bi,bj) = 0. _d 0
115     #endif
116     #ifdef ALLOW_AUTODIFF_TAMC
117     iceFrac(i,j) = 0.
118     #endif
119     ENDDO
120     ENDDO
121    
122     ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
123     snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
124     DO j = jMin, jMax
125     DO i = iMin, iMax
126     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
127     C-- Snow aging :
128     snowAge(i,j,bi,bj) = thSIce_deltaT
129     & + snowAge(i,j,bi,bj)*ageFac
130     IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
131     & snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
132     & * EXP( - snowFac*snowPrc(i,j,bi,bj) )
133     c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
134     c & /hNewSnowAge )
135     C-------
136     C note: Any flux of mass (here fresh water) that enter or leave the system
137     C with a non zero energy HAS TO be counted: add snow precip.
138     icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
139     & - Lfresh*snowPrc(i,j,bi,bj)
140     C--
141     ENDIF
142     ENDDO
143     ENDDO
144    
145     #ifdef ALLOW_DIAGNOSTICS
146     IF ( useDiagnostics ) THEN
147     tmpFac = 1. _d 0
148     CALL DIAGNOSTICS_FRACT_FILL(
149     I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
150     I 0,1,1,bi,bj,myThid)
151     CALL DIAGNOSTICS_FRACT_FILL(
152     I siceAlb, iceMask,tmpFac,1,'SIalbedo',
153     I 0,1,1,bi,bj,myThid)
154     ENDIF
155     #endif /* ALLOW_DIAGNOSTICS */
156     #ifndef ALLOW_ATM2D
157     DO j = jMin, jMax
158     DO i = iMin, iMax
159     siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
160     ENDDO
161     ENDDO
162     #endif
163    
164     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
165     C part.2 : ice-covered fraction ;
166     C change in ice/snow thickness and ice-fraction
167     C note: can only reduce the ice-fraction but not increase it.
168     C-------
169     DO j = jMin, jMax
170     DO i = iMin, iMax
171    
172     tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
173     cphm = cpwater*rhosw*hOceMxL(i,j,bi,bj)
174     frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
175     & * cphm/ocean_deltaT
176     iceFrac(i,j) = iceMask(i,j,bi,bj)
177     flx2oc(i,j) = icFlxSW(i,j,bi,bj)
178     if (frzmltMxL(i,j).GT.0.D0) THEN
179     if ((i.eq.4).and.(j.eq.43)) THEN
180     PRINT *,'in frzmlt loop at:',i,j
181     print *, frzmltMxL(i,j),iceFrac(i,j),
182     & flx2oc(i,j),tOceMxL(i,j,bi,bj)
183     endif
184     endif
185     C-------
186     #ifdef ALLOW_DBUG_THSICE
187     IF ( dBug(i,j,bi,bj) ) THEN
188     IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
189     WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
190     WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
191     & iceFrac(i,j), iceHeight(i,j,bi,bj),
192     & snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
193     WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
194     & tOceMxL(i,j,bi,bj), tFrzOce(i,j),
195     & frzmltMxL(i,j), Qnet(i,j,bi,bj)
196     ENDIF
197     IF (iceFrac(i,j).GT.0.)
198     & WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
199     & iceFrac(i,j), icFlxAtm(i,j,bi,bj),
200     & icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
201     ENDIF
202     #endif
203     ENDDO
204     ENDDO
205    
206     CALL THSICE_CALC_THICKN(
207     I bi, bj, siLo, siHi, sjLo, sjHi,
208     I iMin,iMax, jMin,jMax, dBugFlag,
209     I iceMask(siLo,sjLo,bi,bj), tFrzOce,
210     I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
211     I snowPrc(siLo,sjLo,bi,bj), prcAtm,
212     I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
213     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
214     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
215     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
216     U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
217     O frw2oc, fsalt,
218     I myTime, myIter, myThid )
219    
220     C-- Net fluxes :
221     DO j = jMin, jMax
222     DO i = iMin, iMax
223     IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
224     C- weighted average net fluxes:
225     icFrac = iceMask(i,j,bi,bj)
226     opFrac= 1. _d 0-icFrac
227     #ifdef ALLOW_ATM2D
228     pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
229     pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
230     sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
231     #else
232     icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
233     & - opFrac*Qnet(i,j,bi,bj)
234     icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
235     & + opFrac*rhofw*EmPmR(i,j,bi,bj)
236     Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
237     EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
238     & + opFrac*EmPmR(i,j,bi,bj)
239     saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
240     #endif
241    
242     #ifdef ALLOW_DBUG_THSICE
243     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
244     & 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
245     & iceFrac(i,j), iceHeight(i,j,bi,bj),
246     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
247     #endif
248    
249     #ifndef ALLOW_ATM2D
250     ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
251     icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
252     icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
253     ELSE
254     icFlxAtm(i,j,bi,bj) = 0. _d 0
255     icFrwAtm(i,j,bi,bj) = 0. _d 0
256     #endif
257     ENDIF
258     ENDDO
259     ENDDO
260    
261     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
262     C part.3 : freezing of sea-water
263     C over ice-free fraction and what is left from ice-covered fraction
264     C-------
265     CALL THSICE_EXTEND(
266     I bi, bj, siLo, siHi, sjLo, sjHi,
267     I iMin,iMax, jMin,jMax, dBugFlag,
268     I frzmltMxL, tFrzOce,
269     I tOceMxL(siLo,sjLo,bi,bj),
270     U iceFrac, iceHeight(siLo,sjLo,bi,bj),
271     U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
272     U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
273     U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
274     O flx2oc, frw2oc, fsalt,
275     I myTime, myIter, myThid )
276    
277     DO j = jMin, jMax
278     DO i = iMin, iMax
279     IF (frzmltMxL(i,j).GT.0. _d 0) THEN
280     C-- Net fluxes :
281     #ifdef ALLOW_ATM2D
282     pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
283     pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
284     sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
285     #else
286     Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
287     EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
288     saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
289     #endif
290    
291     #ifdef ALLOW_DBUG_THSICE
292     IF (dBug(i,j,bi,bj)) WRITE(6,1010)
293     & 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
294     & iceFrac(i,j), iceHeight(i,j,bi,bj),
295     & snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
296     #endif
297     ENDIF
298    
299     IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
300     & isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
301     & .AND. iceFrac(i,j) .LE.0. _d 0
302     IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
303     iceMask(i,j,bi,bj)=iceFrac(i,j)
304     IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
305     & snowAge(i,j,bi,bj) = 0. _d 0
306     ELSE
307     iceMask(i,j,bi,bj) = 0. _d 0
308     iceHeight(i,j,bi,bj)= 0. _d 0
309     snowHeight(i,j,bi,bj)=0. _d 0
310     snowAge(i,j,bi,bj) = 0. _d 0
311     Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
312     Tice1(i,j,bi,bj) = 0. _d 0
313     Tice2(i,j,bi,bj) = 0. _d 0
314     Qice1(i,j,bi,bj) = 0. _d 0
315     Qice2(i,j,bi,bj) = 0. _d 0
316     ENDIF
317    
318     #ifdef ATMOSPHERIC_LOADING
319     C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
320     #ifndef ALLOW_ATM2D
321     sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
322     & + iceHeight(i,j,bi,bj)*rhoi
323     & )*iceMask(i,j,bi,bj)
324     #endif
325     #endif
326    
327     ENDDO
328     ENDDO
329    
330     #ifdef ALLOW_BULK_FORCE
331     IF ( useBulkForce ) THEN
332     CALL BULKF_FLUX_ADJUST(
333     I bi, bj, iMin, iMax, jMin, jMax,
334     I isIceFree, myTime, myIter, myThid )
335     ENDIF
336     #endif /* ALLOW_BULK_FORCE */
337    
338     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
339     #endif /* ALLOW_THSICE */
340    
341     RETURN
342     END

  ViewVC Help
Powered by ViewVC 1.1.22