forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
thsice_extend.F
257 lines (238 loc) · 9.1 KB
/
thsice_extend.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
#include "THSICE_OPTIONS.h"
CBOP
C !ROUTINE: THSICE_EXTEND
C !INTERFACE:
SUBROUTINE THSICE_EXTEND(
I bi, bj,
I iMin,iMax, jMin,jMax, dBugFlag,
I fzMlOc, tFrz, tOce,
U icFrac, hIce, hSnow,
U tSrf, tIc1, tIc2, qIc1, qIc2,
O flx2oc, frw2oc, fsalt,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R THSICE_EXTEND
C | o Extend sea-ice area incresing ice fraction
C *==========================================================*
C | o incorporate surplus of energy to
C | make new ice or make ice grow laterally
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C bi,bj :: tile indices
C iMin,iMax :: computation domain: 1rst index range
C jMin,jMax :: computation domain: 2nd index range
C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
C--- Input:
C iceMask :: sea-ice fractional mask [0-1]
C fzMlOc (esurp) :: ocean mixed-layer freezing/melting potential [W/m2]
C tFrz (Tf) :: sea-water freezing temperature [oC] (function of S)
C tOce :: surface level oceanic temperature [oC]
C--- Modified (input&output):
C icFrac(iceFrac):: fraction of grid area covered in ice
C hIce (iceThick):: ice height [m]
C hSnow :: snow height [m]
C tSrf :: surface (ice or snow) temperature [oC]
C tIc1 :: temperature of ice layer 1 [oC]
C tIc2 :: temperature of ice layer 2 [oC]
C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level
C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level
C--- Output
C flx2oc (=) :: (additional) heat flux to ocean [W/m2] (+=dwn)
C frw2oc (=) :: (additional) fresh water flux to ocean [kg/m2/s] (+=dwn)
C fsalt (=) :: (additional) salt flux to ocean [g/m2/s] (+=dwn)
C--- Input:
C myTime :: current Time of simulation [s]
C myIter :: current Iteration number in simulation
C myThid :: my Thread Id number
INTEGER bi,bj
INTEGER iMin, iMax
INTEGER jMin, jMax
LOGICAL dBugFlag
c _RL iceMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL fzMlOc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tFrz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL icFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hSnow (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tSrf (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tIc1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tIc2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL qIc1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL qIc2 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_THSICE
C !LOCAL VARIABLES:
C--- local copy of input/output argument list variables (see description above)
_RL esurp
_RL Tf
_RL iceFrac
_RL iceThick
_RL qicen(nlyr)
C == Local variables ==
C iceVol :: previous ice volume
C newIce :: new ice volume to produce
C hNewIce :: thickness of new ice to form
C iceFormed :: ice-volume formed (new ice volume = iceVol+iceFormed )
C qicAv :: mean enthalpy of ice (layer 1 & 2) [J/m^3]
_RL deltaTice ! time-step for ice model
_RL iceVol
_RL newIce
_RL hNewIce
_RL iceFormed
_RL qicAv
INTEGER i,j ! loop indices
C- define grid-point location where to print debugging values
#include "THSICE_DEBUG.h"
1020 FORMAT(A,1P4E11.3)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
deltaTice = thSIce_deltaT
#ifdef ALLOW_AUTODIFF_TAMC
DO j = 1-OLy, sNy+OLy
DO i = 1-OLx, sNx+OLx
flx2oc(i,j) = 0. _d 0
frw2oc(i,j) = 0. _d 0
fsalt (i,j) = 0. _d 0
ENDDO
ENDDO
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ticekey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hIce(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE hSnow(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE icFrac(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE qIc1(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE qIc2(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
#endif
DO j = jMin, jMax
DO i = iMin, iMax
c#ifdef ALLOW_AUTODIFF_TAMC
c ikey_1 = i
c & + sNx*(j-1)
c & + sNx*sNy*act1
c & + sNx*sNy*max1*act2
c & + sNx*sNy*max1*max2*act3
c & + sNx*sNy*max1*max2*max3*act4
c#endif /* ALLOW_AUTODIFF_TAMC */
c#ifdef ALLOW_AUTODIFF_TAMC
cCADJ STORE hice(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE hsnow(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE icfrac(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE qic1(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE qic2(i,j) = comlev1_thsice_1, key=ikey_1
c#endif
IF (fzMlOc(i,j).GT.0. _d 0) THEN
esurp = fzMlOc(i,j)
Tf = tFrz(i,j)
iceFrac = icFrac(i,j)
iceThick= hIce(i,j)
qicen(1)= qIc1(i,j)
qicen(2)= qIc2(i,j)
C---
C-- start ice
iceFormed = 0. _d 0
iceVol = iceFrac*iceThick
C- enthalpy of new ice to form :
IF ( iceFrac.LE.0. _d 0 ) THEN
qicen(1)= -cpWater*Tmlt1
& + cpIce *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)
qicen(2)= -cpIce *Tf + Lfresh
ENDIF
qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0
newIce = esurp*deltaTice/qicAv
IF ( icFrac(i,j).EQ.0. _d 0 ) THEN
C- to keep identical results (as it use to be):
c-old_v IF ( newIce.GE.hThinIce*iceMaskMin ) THEN
C- here we allow to form ice earlier (as soon as min-ice-vol is reached)
c-new_v:
IF ( newIce.GT.hIceMin*iceMaskMin ) THEN
C- if there is no ice in grid-cell and enough ice to form:
C- make ice over iceMaskMin fraction, up to hThinIce,
C and if more ice to form, then increase fraction
iceThick = MIN(hThinIce,newIce/iceMaskMin)
iceThick = MAX(iceThick,newIce/iceMaskMax)
iceFrac = newIce/iceThick
iceFormed = newIce
ENDIF
ELSEIF ( iceVol.LT.hiMax*iceMaskMax ) THEN
C- if there is already some ice
C create ice with same thickness or hNewIceMax (the smallest of the 2)
hNewIce = MIN(iceThick,hNewIceMax)
iceFrac = MIN(icFrac(i,j)+newIce/hNewIce,iceMaskMax)
C- update thickness: area weighted average
c-new_v:
iceThick = MIN(hiMax,(iceVol+newIce)/iceFrac)
C- to keep identical results: comment the line above and uncomment line below:
c-old_v iceFrac = MIN(icFrac(i,j)+newIce/iceThick,iceMaskMax)
iceFormed = iceThick*iceFrac - iceVol
C- spread snow out over ice
hSnow(i,j) = hSnow(i,j)*icFrac(i,j)/iceFrac
ENDIF
C- oceanic fluxes:
flx2oc(i,j)= qicAv*iceFormed/deltaTice
frw2oc(i,j)= -rhoi*iceFormed/deltaTice
fsalt(i,j)= -(rhoi*saltIce)*iceFormed/deltaTice
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) THEN
WRITE(6,1020) 'ThSI_EXT: iceH, newIce, newIceFrac=',
& iceThick, newIce, iceFrac-icFrac(i,j)
WRITE(6,1020) 'ThSI_EXT: iceFrac,flx2oc,fsalt,frw2oc=',
& iceFrac,flx2oc(i,j),fsalt(i,j),frw2oc(i,j)
ENDIF
#endif
#ifdef CHECK_ENERGY_CONSERV
CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 1,
I icFrac(i,j), iceFrac, iceThick, hSnow(i,j), qicen,
I flx2oc(i,j), frw2oc(i,j), fsalt(i,j),
I myTime, myIter, myThid )
#endif /* CHECK_ENERGY_CONSERV */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Update Sea-Ice state output:
IF ( iceFrac.GT.0. _d 0 .AND. icFrac(i,j).EQ.0. _d 0) THEN
c hSnow(i,j) = 0. _d 0
tSrf(i,j) = tFrz(i,j)
tIc1(i,j) = tFrz(i,j)
tIc2(i,j) = tFrz(i,j)
qIc1(i,j) = qicen(1)
qIc2(i,j) = qicen(2)
ENDIF
icFrac(i,j) = iceFrac
hIce(i,j) = iceThick
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_THSICE */
RETURN
END