forked from NOAA-EMC/HYCOM-src
-
Notifications
You must be signed in to change notification settings - Fork 0
/
icloan.F90
395 lines (394 loc) · 16.6 KB
/
icloan.F90
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
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
#if defined (USE_NUOPC_CESMBETA) || (ESPC_COUPLE)
#define USE_NUOPC_GENERIC 1
#endif
#if defined(ROW_LAND)
#define SEA_P .true.
#define SEA_U .true.
#define SEA_V .true.
#elif defined(ROW_ALLSEA)
#define SEA_P allip(j).or.ip(i,j).ne.0
#define SEA_U alliu(j).or.iu(i,j).ne.0
#define SEA_V alliv(j).or.iv(i,j).ne.0
#else
#define SEA_P ip(i,j).ne.0
#define SEA_U iu(i,j).ne.0
#define SEA_V iv(i,j).ne.0
#endif
subroutine icloan(m,n)
use mod_xc ! HYCOM communication interface
use mod_cb_arrays ! HYCOM saved arrays
#if defined (USE_NUOPC_GENERIC)
use mod_import ! HYCOM merge import
#endif
implicit none
!
integer m,n
!
! --- 'energy loan' ice model. no advection, no dynamics. ice amount
! --- represents energy 'loaned' to water column to prevent wintertime
! --- cooling below freezing level. loan is paid back in summer.
!
! --- modified version for ice-ocean "coupling".
! --- freeze/melt energy from relaxation to the freezing temperature.
!
integer i,j
real tfrz,tsur,tmxl,smxl,hfrz,paybak,borrow,hice,thkimx,t2f
real radfl,swfl,tdif,wind,airt,pair,rair,snsibl,emnp,dtrmui
real thkimxy(jdm)
!
! --- hice = actual ice thickness (m), local variable
!
! --- thkice = average ice thickness, i.e. hice x covice (m)
! --- covice = ice coverage, i.e. cell fraction (0.0 to 1.0)
! --- temice = ice surface temperature (degC)
! --- flxice = cell average heat flux under ice (W/m^2) into ocean
! --- fswice = cell average swv flux under ice (W/m^2) into ocean
! --- sflice = cell average salt flux under ice (psu kg/m^2/s) into ocean
! --- wflice = cell average water flux under ice ( kg/m^2/s) into ocean
! --- wflfrz = cell average water flux under ice due to freezing/melting
!
! --- icefrq = e-folding time scale back to tfrz (no. time steps)
! --- thkfrz = maximum thickness of near-surface freezing zone (m)
! --- tfrz_0 = ice melting point (degC) at S=0psu
! --- tfrz_s = gradient of ice melting point (degC/psu)
! --- ticegr = vertical temperature gradient inside ice (deg/m)
! --- (0.0 to get ice surface temp. from atmos. surtmp)
! --- hicemn = minimum ice thickness (m)
! --- hicemx = maximum ice thickness (m)
!
real tfrz_n,ticemn,ticemx,salice,rhoice,fusion,meltmx
parameter (tfrz_n= -1.79, & ! nominal ice melting point (degC)
ticemn=-50.0, & ! minimum ice surface temperature (degC)
ticemx= 0.0, & ! maximum ice surface temperature (degC)
salice= 4.0, & ! salinity of ice (psu) - same as CICE
rhoice=917.0, & ! density of ice (kg/m**3)
fusion=334.e3, & ! latent heat of fusion (J/kg)
meltmx= 33.e-7)! max. ice melting rate (m/s), 0.285 m/day
!
real fluxmx !max. ice melting flux (W/m^2)
parameter (fluxmx=meltmx*fusion*rhoice) !~1000 W/m^2 - like CICE
!
real csice,csubp,pairc,rgas,tzero
parameter (csice =0.0006, & !ice-air sensible exchange coefficient
csubp =1005.7, & !specific heat of air (j/kg/deg)
pairc=1013.0*100.0, & !air pressure (mb) * 100
rgas =287.1, & !gas constant (j/kg/k)
tzero=273.16) !celsius to kelvin temperature offset
!
real sb_cst
parameter (sb_cst=5.67e-8) !Stefan-Boltzman constant, for lwflag=-1
!
!# include "stmt_fns.h"
!
dtrmui = baclin/(1.0*86400.0) !dt*1/1days
!
! --- energy loan: add extra energy to the ocean to keep SST from dropping
! --- below tfrz in winter. return this borrowed energy to the 'energy bank'
! --- in summer.
!
! --- salt loan: analogous to energy loan.
!
!$OMP PARALLEL DO PRIVATE(j,i, &
!$OMP t2f,hfrz,smxl,tmxl,tfrz,borrow,paybak) &
!$OMP SCHEDULE(STATIC,jblk)
do j=1,jj
thkimxy(j)=0.0 !simplifies OpenMP parallelization
do i=1,ii
if (SEA_P) then
if (ishlf(i,j).eq.0) then !under an ice shelf
flxice(i,j)=0.0
wflice(i,j)=0.0
sflice(i,j)=0.0
thkice(i,j)=hicemx
util1(i,j)=0.0
else !standard ocean point
! --- relax to tfrz with e-folding time of icefrq time steps
! --- assuming the effective surface layer thickness (hfrz)
! --- is at most thkfrz meters
! --- multiply by dpbl(i,j)/hfrz to get the actual e-folding time
! --- icefrq==1 is "Instant Relaxation" and an e-folding time
! --- of 30 days is consistent with the "Drag Law" approach
! --- D.M. Holland (1998) On the Parameterization of Basal
! --- Heat Flux for Sea-ice Modelling; Geophysica 34 pp 1-21
! --- when coupling, icefrq must be at least the coupling interval
!
hfrz = min( thkfrz*onem, dpbl(i,j) )
t2f = (spcifh*hfrz)/(baclin*icefrq*g)
smxl = saln(i,j,1,n)
tmxl = temp(i,j,1,n)
tfrz = tfrz_0 + smxl*tfrz_s !salinity dependent freezing point
borrow = (tfrz-tmxl)*t2f !W/m^2 into ocean
!
! --- limit heat flux range (for both forming and melting ice)
borrow=max( -fluxmx, min( fluxmx, borrow ) )
!
!diag if (i.eq.itest .and. j.eq.jtest) then
!diag write (lp,'(i9,2i5,a,5f9.3)') &
!diag nstep,i+i0,j+j0,' t,tfrz,flx,hfrz,cov:', &
!diag tmxl,tfrz,borrow,hfrz*qonem,covice(i,j)
!diag endif
!
if (tmxl.lt.tfrz) then
!
! --- add energy to move tmxl towards tfrz (only if tmxl < tfrz)
! --- include some dependance on sea ice coverage
!
flxice(i,j)=borrow*max(covice(i,j),0.1) !+ve
wflice(i,j)= -flxice(i,j)/fusion !-ve
thkice(i,j)=thkice(i,j)-wflice(i,j)*(baclin/rhoice) !+ve inc.
! --- brine rejection as ice forms
sflice(i,j)= wflice(i,j)*min(smxl,salice) !-ve
elseif (thkice(i,j).gt.0.0) then !tmxl > tfrz
!
! --- ice, so return the borrowed amount whenever tmxl > tfrz
! --- only over the fraction of the cell where ice exists
!
paybak=min( -borrow*covice(i,j), & !+ve
thkice(i,j)*(fusion*rhoice/baclin) )
flxice(i,j)= -paybak !-ve
wflice(i,j)= -flxice(i,j)/fusion !+ve
thkice(i,j)=thkice(i,j)-wflice(i,j)*(baclin/rhoice) !-ve inc.
! --- brine recovery from melting ice
sflice(i,j)= wflice(i,j)*salice !+ve
else !tmxl > tfrz & thkice(i,j) == 0.0
!
! --- no ice.
!
flxice(i,j)=0.0
wflice(i,j)=0.0
sflice(i,j)=0.0
!
if (icmflg.eq.2) then
!
! --- add extra cooling under the ice mask (tsur<=tfrz_n)
! --- don't allow a new tsur maximum, to preserve sea ice
!
if (natm.eq.2) then
tsur = min( max( surtmp(i,j,l0), surtmp(i,j,l1) ), &
surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1 )
elseif (yrflag.lt.2) then
tsur = min( max( surtmp(i,j,l0), surtmp(i,j,l1), &
surtmp(i,j,l2), surtmp(i,j,l3) ), &
surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1+ &
surtmp(i,j,l2)*w2+surtmp(i,j,l3)*w3 )
else
tsur = min( max( surtmp(i,j,l0), surtmp(i,j,l1) ), &
surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1 )
endif
if (tsur.le.tfrz_n) then
surflx(i,j)=surflx(i,j)+borrow*covice(i,j)
endif
endif !icmflg.eq.2
endif
!
util1(i,j)=max(thkice(i,j)-hicemx,0.0) !icex = ice exceeding hicemx
thkimxy(j)=max(thkimxy(j),thkice(i,j))
endif !ishlf:else
wflfrz(i,j)=wflice(i,j) !diagnostic only
endif !ip
enddo !i
enddo !j
!$OMP END PARALLEL DO
!
thkimx=maxval(thkimxy(1:jj))
call xcmaxr(thkimx)
!
! --- spread out portion of ice thicker than hicemx
if (thkimx.gt.hicemx) then
call psmooth(util1, 0,0, ishlf, util2) !smooth icex
endif
!
!$OMP PARALLEL DO PRIVATE(j,i,hice,smxl,tfrz, &
!$OMP radfl,tdif,wind,airt,rair,snsibl,emnp) &
!$OMP SCHEDULE(STATIC,jblk)
do j=1,jj
do i=1,ii
if (SEA_P) then
if (ishlf(i,j).eq.0) then !under an ice shelf
thkice(i,j)=hicemx
covice(i,j)=1.0
fswice(i,j)=0.0
temice(i,j)=ticemn
else !standard ocean point
thkice(i,j)=util1(i,j)+min(thkice(i,j),hicemx) !icex_sm+rest
!
! --- compute fractional ice coverage for energy flux calculation
if (thkice(i,j).lt.1.e-5*hicemn) then
covice(i,j)=0.0
else
covice(i,j)=min(1.0,thkice(i,j)*(1.0/hicemn))
hice=thkice(i,j)/covice(i,j) !minimum of hicemn
end if
!
if (icmflg.eq.3) then
! --- relax to sea ice concentration from coupler
! --- ice thickness is therefore always then 0 and hicemn
covice(i,j)=covice(i,j)+dtrmui*(si_c(i,j)-covice(i,j))
thkice(i,j)=covice(i,j)*hicemn
hice=hicemn
endif
!
! --- compute ice surface temperature
if (covice(i,j).eq.0.0) then
temice(i,j)=ticemx
elseif (icmflg.eq.3 .and. si_c(i,j).gt.0.0) then !from coupler
temice(i,j)=max( ticemn, min( ticemx, si_t(i,j) ) )
elseif (ticegr.eq.0.0) then !use surtmp
if (natm.eq.2) then
temice(i,j)=max( ticemn, &
min( ticemx, &
surtmp(i,j,l0)*w0+ &
surtmp(i,j,l1)*w1 ) )
else
temice(i,j)=max( ticemn, &
min( ticemx, &
surtmp(i,j,l0)*w0+ &
surtmp(i,j,l1)*w1+ &
surtmp(i,j,l2)*w2+ &
surtmp(i,j,l3)*w3 ) )
endif !natm
else
temice(i,j)=max( ticemn, ticemx-ticegr*hice )
endif
!
! --- atmosphere to ice surface exchange is applied to the ocean,
! --- i.e. use the "energy-loan" approach.
! --- don't apply to the ocean when coupling, because it has
! --- already been applied to the ice.
!
if (icmflg.ne.3 .and. covice(i,j).gt.0.0) then
! --- net radiative thermal flux (w/m**2) +ve into ocean/ice
if (lwflag.gt.-1) then
! --- radflx's Qsw includes the atmos. model's surface albedo,
! --- i.e. it already allows for ice&snow where it is observed.
if (natm.eq.2) then
radfl=radflx(i,j,l0)*w0+radflx(i,j,l1)*w1
else
radfl=radflx(i,j,l0)*w0+radflx(i,j,l1)*w1 &
+radflx(i,j,l2)*w2+radflx(i,j,l3)*w3
endif !natm
if (lwflag.gt.1) then
! --- longwave correction to radfl (Qsw+Qlw).
! --- this will be ~zero for ticegr==0.0 (temice=surtmp)
if (natm.eq.2) then
tdif = temice(i,j) - &
( surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1 )
else
tdif = temice(i,j) - &
( surtmp(i,j,l0)*w0+surtmp(i,j,l1)*w1 &
+surtmp(i,j,l2)*w2+surtmp(i,j,l3)*w3 )
endif !natm
!correction is blackbody radiation from tdif at temice
radfl = radfl - (4.506+0.0554*temice(i,j)) * tdif
endif
else !lwflag.eq.-1
! --- input radflx is Qlwdn
! --- input swflx is net Qsw with ocean albedo (~ 1-0.09)
! --- convert swfl to net Qsw with sea ice albedo (1-0.6)
! --- i.e. scale by(1-0.6)/(1-0.09)=0.44
! --- convert radfl to net Qlw + Qsw where sea ice
if (natm.eq.2) then
radfl=radflx(i,j,l0)*w0+radflx(i,j,l1)*w1
swfl=swflx (i,j,l0)*w0+swflx (i,j,l1)*w1
else
radfl=radflx(i,j,l0)*w0+radflx(i,j,l1)*w1 &
+radflx(i,j,l2)*w2+radflx(i,j,l3)*w3
swfl=swflx (i,j,l0)*w0+swflx (i,j,l1)*w1 &
+swflx (i,j,l2)*w2+swflx (i,j,l3)*w3
endif !natm
radfl = radfl - sb_cst*(temice(i,j)+tzero)**4 &
+ 0.44d0*swfl !!Alex
endif !lwflux
!
if (flxflg.ne.3) then
if (natm.eq.2) then
! --- wind speed (m/s)
wind=wndspd(i,j,l0)*w0+wndspd(i,j,l1)*w1
! --- air temperature (C)
airt=airtmp(i,j,l0)*w0+airtmp(i,j,l1)*w1
else
! --- wind speed (m/s)
wind=wndspd(i,j,l0)*w0+wndspd(i,j,l1)*w1 &
+wndspd(i,j,l2)*w2+wndspd(i,j,l3)*w3
! --- air temperature (C)
airt=airtmp(i,j,l0)*w0+airtmp(i,j,l1)*w1 &
+airtmp(i,j,l2)*w2+airtmp(i,j,l3)*w3
endif !natm
if (mslprf .or. flxflg.eq.6) then
if (natm.eq.2) then
pair=mslprs(i,j,l0)*w0+mslprs(i,j,l1)*w1 &
+prsbas
else
pair=mslprs(i,j,l0)*w0+mslprs(i,j,l1)*w1 &
+mslprs(i,j,l2)*w2+mslprs(i,j,l3)*w3 &
+prsbas
endif !natm
else
pair=pairc
endif
rair = pair/(rgas*(tzero+airt))
snsibl = csubp*rair*wind*csice*(temice(i,j)-airt)
#if defined (USE_NUOPC_GENERIC)
if (cpl_merge) snsibl=hycom_imp_mrg_sensflx(i,j,snsibl)
#endif
else
snsibl = 0.0 !already in total flux (i.e. in radfl)
endif
flxice(i,j) = flxice(i,j) + &
covice(i,j)*(radfl - snsibl) !no evap
!
! --- add a time-invarient net heat flux offset
if (flxoff) then
flxice(i,j) = flxice(i,j) + covice(i,j)*offlux(i,j)
endif
!
! --- emnp = evaporation minus precipitation (m/sec) into atmos.
! --- no evap (sublimation) over ice, all precip enters ocean
if (pcipf) then
if (natm.eq.2) then
emnp = -( precip(i,j,l0)*w0+precip(i,j,l1)*w1)
else
emnp = -( precip(i,j,l0)*w0+precip(i,j,l1)*w1 &
+precip(i,j,l2)*w2+precip(i,j,l3)*w3)
endif !natm
else
emnp = 0.0
endif
! --- wflice = water flux (m/s kg/m**2/sec) into ocean under ice
wflice(i,j) = wflice(i,j) - covice(i,j)*emnp*rhoref
endif !covice
fswice(i,j) = 0.0 !no penetrating Qsw under ice
endif !ishlf:else
endif !ip
enddo !i
enddo !j
!$OMP END PARALLEL DO
!
return
end subroutine icloan
!
!
!> Revision history
!>
!> June 2000 - conversion to SI units
!> July 2000 - switched sign convention for vertical fluxes (now >0 if down)
!> May 2003 - added option to impose an ice mask
!> June 2003 - added 8 time step e-folding time scale
!> June 2003 - limited rate of ice formation
!> June 2003 - replaced constant saldif with smxl-salice
!> Mar. 2005 - freezing point linearly dependent on salinity
!> Mar. 2005 - ice surface temperature optionally from surtmp
!> June 2006 - modified version for ice-ocean "coupling"
!> Nov. 2011 - don't apply atmosphere to ice surface exchange when "coupling"
!> May 2012 - limit brine rejection to be a non-negative salt flux
!> July 2012 - flxice and sflice now correctly represent cell average under ice
!> Nov. 2012 - weaker dependance on covice when freezing
!> Jan. 2014 - added natm
!> Apr. 2014 - added pair for time varying msl pressure (mslprf)
!> Apr. 2014 - added ice shelf logic (ishlf)
!> May 2014 - use land/sea masks (e.g. ip) to skip land
!> Aug. 2018 - added sflfrz (now wflfrz) as a diagnostic field
!> Aug. 2018 - icloan is not leapfrog (baclin, not delt1)
!> Nov. 2018 - virtual salt flux replaced with water and actual salt flux
!> Nov. 2018 - added lwflag=-1 for input radflx=Qlwdn
!> Nov. 2018 - allow for difference in ocean and sea ice albedo when lwflag=-1