/[MITgcm]/MITgcm/pkg/thsice/thsice_extend.F
ViewVC logotype

Contents of /MITgcm/pkg/thsice/thsice_extend.F

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


Revision 1.5 - (show annotations) (download)
Mon Apr 16 22:38:24 2007 UTC (17 years, 2 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint59
Changes since 1.4: +34 -2 lines
First set of modifs for TAF-ing thsice.

1 C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_extend.F,v 1.4 2006/05/25 18:03:24 jmc Exp $
2 C $Name: $
3
4 #include "THSICE_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: THSICE_EXTEND
8 C !INTERFACE:
9 SUBROUTINE THSICE_EXTEND(
10 I bi, bj, siLo, siHi, sjLo, sjHi,
11 I iMin,iMax, jMin,jMax, dBugFlag,
12 I fzMlOc, tFrz, tOce,
13 U icFrac, hIce, hSnow,
14 U tSrf, tIc1, tIc2, qIc1, qIc2,
15 O flx2oc, frw2oc, fsalt,
16 I myTime, myIter, myThid )
17 C !DESCRIPTION: \bv
18 C *==========================================================*
19 C | S/R THSICE_EXTEND
20 C | o Extend sea-ice area incresing ice fraction
21 C *==========================================================*
22 C | o incorporate surplus of energy to
23 C | make new ice or make ice grow laterally
24 C *==========================================================*
25 C \ev
26
27 C !USES:
28 IMPLICIT NONE
29
30 C == Global variables ==
31 #include "EEPARAMS.h"
32 #include "THSICE_SIZE.h"
33 #include "THSICE_PARAMS.h"
34 #ifdef ALLOW_AUTODIFF_TAMC
35 # include "SIZE.h"
36 # include "tamc.h"
37 # include "tamc_keys.h"
38 #endif
39
40 C !INPUT/OUTPUT PARAMETERS:
41 C == Routine Arguments ==
42 C siLo,siHi :: size of input/output array: 1rst dim. lower,higher bounds
43 C sjLo,sjHi :: size of input/output array: 2nd dim. lower,higher bounds
44 C bi,bj :: tile indices
45 C iMin,iMax :: computation domain: 1rst index range
46 C jMin,jMax :: computation domain: 2nd index range
47 C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
48 C--- Input:
49 C iceMask :: sea-ice fractional mask [0-1]
50 C fzMlOc (esurp) :: ocean mixed-layer freezing/melting potential [W/m2]
51 C tFrz (Tf) :: sea-water freezing temperature [oC] (function of S)
52 C tOce (sst) :: surface level oceanic temperature [oC]
53 C--- Modified (input&output):
54 C icFrac(iceFrac):: fraction of grid area covered in ice
55 C hIce (iceThick):: ice height [m]
56 C hSnow(snowThick):: snow height [m]
57 C tSrf :: surface (ice or snow) temperature [oC]
58 C tIc1 :: temperature of ice layer 1 [oC]
59 C tIc2 :: temperature of ice layer 2 [oC]
60 C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level
61 C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level
62 C--- Output
63 C flx2oc (=) :: (additional) heat flux to ocean [W/m2] (+=dwn)
64 C frw2oc (=) :: (additional) fresh water flux to ocean [kg/m2/s] (+=dwn)
65 C fsalt (=) :: (additional) salt flux to ocean [g/m2/s] (+=dwn)
66 C--- Input:
67 C myTime :: current Time of simulation [s]
68 C myIter :: current Iteration number in simulation
69 C myThid :: my Thread Id number
70 INTEGER siLo, siHi, sjLo, sjHi
71 INTEGER bi,bj
72 INTEGER iMin, iMax
73 INTEGER jMin, jMax
74 LOGICAL dBugFlag
75 c _RL iceMask(siLo:siHi,sjLo:sjHi)
76 _RL fzMlOc (siLo:siHi,sjLo:sjHi)
77 _RL tFrz (siLo:siHi,sjLo:sjHi)
78 _RL tOce (siLo:siHi,sjLo:sjHi)
79 _RL icFrac (siLo:siHi,sjLo:sjHi)
80 _RL hIce (siLo:siHi,sjLo:sjHi)
81 _RL hSnow (siLo:siHi,sjLo:sjHi)
82 _RL tSrf (siLo:siHi,sjLo:sjHi)
83 _RL tIc1 (siLo:siHi,sjLo:sjHi)
84 _RL tIc2 (siLo:siHi,sjLo:sjHi)
85 _RL qIc1 (siLo:siHi,sjLo:sjHi)
86 _RL qIc2 (siLo:siHi,sjLo:sjHi)
87 _RL flx2oc (siLo:siHi,sjLo:sjHi)
88 _RL frw2oc (siLo:siHi,sjLo:sjHi)
89 _RL fsalt (siLo:siHi,sjLo:sjHi)
90 _RL myTime
91 INTEGER myIter
92 INTEGER myThid
93 CEOP
94
95 #ifdef ALLOW_THSICE
96 C !LOCAL VARIABLES:
97 C--- local copy of input/output argument list variables (see description above)
98 _RL esurp
99 _RL Tf
100 _RL sst
101 _RL iceFrac
102 _RL iceThick
103 _RL snowThick
104 _RL qicen(nlyr)
105
106 C == Local variables ==
107 C qicAv :: mean enthalpy of ice (layer 1 & 2) [J/m^3]
108 _RL deltaTice ! time-step for ice model
109 _RL newIce
110 _RL newIceFrac
111 _RL qicAv
112 INTEGER i,j ! loop indices
113
114 C- define grid-point location where to print debugging values
115 #include "THSICE_DEBUG.h"
116
117 1010 FORMAT(A,I3,3F8.3)
118 1020 FORMAT(A,1P4E11.3)
119
120 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
121
122 #ifdef ALLOW_AUTODIFF_TAMC
123 act1 = bi - myBxLo(myThid)
124 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
125 act2 = bj - myByLo(myThid)
126 max2 = myByHi(myThid) - myByLo(myThid) + 1
127 act3 = myThid - 1
128 max3 = nTx*nTy
129 act4 = ikey_dynamics - 1
130 #endif /* ALLOW_AUTODIFF_TAMC */
131
132 deltaTice = thSIce_deltaT
133
134 DO j = jMin, jMax
135 DO i = iMin, iMax
136 #ifdef ALLOW_AUTODIFF_TAMC
137 ikey_1 = i
138 & + sNx*(j-1)
139 & + sNx*sNy*act1
140 & + sNx*sNy*max1*act2
141 & + sNx*sNy*max1*max2*act3
142 & + sNx*sNy*max1*max2*max3*act4
143 #endif /* ALLOW_AUTODIFF_TAMC */
144 C--
145 #ifdef ALLOW_AUTODIFF_TAMC
146 CADJ STORE hice(i,j) = comlev1_thsice_1, key=ikey_1
147 CADJ STORE hsnow(i,j) = comlev1_thsice_1, key=ikey_1
148 CADJ STORE icfrac(i,j) = comlev1_thsice_1, key=ikey_1
149 CADJ STORE qic1(i,j) = comlev1_thsice_1, key=ikey_1
150 CADJ STORE qic2(i,j) = comlev1_thsice_1, key=ikey_1
151 #endif
152
153 IF (fzMlOc(i,j).GT.0. _d 0) THEN
154 esurp = fzMlOc(i,j)
155 Tf = tFrz(i,j)
156 sst = tOce(i,j)
157 iceFrac = icFrac(i,j)
158 iceThick= hIce(i,j)
159 snowThick=hSnow(i,j)
160 qicen(1)= qIc1(i,j)
161 qicen(2)= qIc2(i,j)
162 C---
163 C-- start ice
164 newIceFrac = 0. _d 0
165
166 C- enthalpy of new ice to form :
167 IF ( iceFrac.LE.0. _d 0 ) THEN
168 qicen(1)= -cpwater*Tmlt1
169 & + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)
170 qicen(2)= -cpice *Tf + Lfresh
171 ENDIF
172 qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0
173 newIce = esurp*deltaTice/qicAv
174
175 IF (icFrac(i,j).EQ.0. _d 0) THEN
176 c IF (newIce.GE.himin*iceMaskmax) THEN
177 C- jmc: above is the original version, but below seems more logical:
178 IF (newIce.GE.himin0*iceMaskmin) THEN
179 C- if there is no ice in grid and enough ice to form:
180 iceThick = MAX(himin0,newIce/iceMaskmax)
181 newIceFrac = MIN(newIce/himin0,iceMaskmax)
182 iceFrac = newIceFrac
183 sst=Tf
184 ENDIF
185 ELSE
186 C- if there is already some ice
187 newIceFrac=MIN(newIce/iceThick,iceMaskmax-icFrac(i,j))
188 iceFrac = icFrac(i,j) + newIceFrac
189 C- spread snow out over ice
190 snowThick = snowThick*icFrac(i,j)/iceFrac
191 sst=(1. _d 0-newIceFrac)*sst+newIceFrac*Tf
192 ENDIF
193 C- oceanic fluxes:
194 flx2oc(i,j)= iceThick*newIceFrac*qicAv/deltaTice
195 frw2oc(i,j)=-(rhoi*iceThick)*newIceFrac/deltaTice
196 fsalt(i,j)= -(rhoi*iceThick*saltice)*newIceFrac/deltaTice
197
198 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
199 #ifdef ALLOW_DBUG_THSICE
200 IF ( dBug(i,j,bi,bj) ) THEN
201 WRITE(6,1020) 'ThSI_EXT: iceH, newIce, newIceFrac=',
202 & iceThick, newIce, newIceFrac
203 WRITE(6,1020) 'ThSI_EXT: iceFrac,flx2oc,fsalt,frw2oc=',
204 & iceFrac,flx2oc(i,j),fsalt(i,j),frw2oc(i,j)
205 ENDIF
206 #endif
207 #ifdef CHECK_ENERGY_CONSERV
208 CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 1,
209 I icFrac(i,j), iceFrac, iceThick, snowThick, qicen,
210 I flx2oc(i,j), frw2oc(i,j), fsalt(i,j),
211 I myTime, myIter, myThid )
212 #endif /* CHECK_ENERGY_CONSERV */
213 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
214 C-- Update Sea-Ice state output:
215 IF ( iceFrac.GT.0. _d 0 .AND. icFrac(i,j).EQ.0. _d 0) THEN
216 tSrf(i,j) = tFrz(i,j)
217 tIc1(i,j) = tFrz(i,j)
218 tIc2(i,j) = tFrz(i,j)
219 qIc1(i,j) = qicen(1)
220 qIc2(i,j) = qicen(2)
221 ENDIF
222 icFrac(i,j) = iceFrac
223 hIce(i,j) = iceThick
224 hSnow(i,j ) = snowThick
225 ENDIF
226 ENDDO
227 ENDDO
228
229 #endif /* ALLOW_THSICE */
230
231 RETURN
232 END

  ViewVC Help
Powered by ViewVC 1.1.22