Newer
Older
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
MODULE icbthm
!!======================================================================
!! *** MODULE icbthm ***
!! Icebergs: thermodynamics routines for icebergs
!!======================================================================
!! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code
!! - ! 2011-03 (Madec) Part conversion to NEMO form
!! - ! Removal of mapping from another grid
!! - ! 2011-04 (Alderson) Split into separate modules
!! - ! 2011-05 (Alderson) Use tmask instead of tmask_i
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! icb_thm : initialise
!! reference for equations - M = Martin + Adcroft, OM 34, 2010
!!----------------------------------------------------------------------
USE par_oce ! NEMO parameters
USE dom_oce ! NEMO domain
USE in_out_manager ! NEMO IO routines, numout in particular
USE lib_mpp ! NEMO MPI routines, ctl_stop in particular
USE phycst ! NEMO physical constants
USE sbc_oce
USE eosbn2 ! equation of state
USE lib_fortran, ONLY : DDPDD
USE icb_oce ! define iceberg arrays
USE icbutl ! iceberg utility routines
USE icbdia ! iceberg budget routines
IMPLICIT NONE
PRIVATE
PUBLIC icb_thm ! routine called in icbstp.F90 module
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: icbthm.F90 15088 2021-07-06 13:03:34Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE icb_thm( kt )
!!----------------------------------------------------------------------
!! *** ROUTINE icb_thm ***
!!
!! ** Purpose : compute the iceberg thermodynamics.
!!
!! ** Method : - See Martin & Adcroft, Ocean Modelling 34, 2010
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! timestep number, just passed to icb_utl_print_berg
!
INTEGER :: ii, ij, jk, ikb
REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb
REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv
REAL(wp) :: zSSS, zfzpt
REAL(wp) :: zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12
REAL(wp) :: zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb
REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw
REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, zzMv
TYPE(iceberg), POINTER :: this, next
TYPE(point) , POINTER :: pt
!
COMPLEX(dp), DIMENSION(jpi,jpj) :: cicb_melt, cicb_hflx
!!----------------------------------------------------------------------
!
!! initialiaze cicb_melt and cicb_heat
cicb_melt = CMPLX( 0.e0, 0.e0, dp )
cicb_hflx = CMPLX( 0.e0, 0.e0, dp )
!
z1_rday = 1._wp / rday
z1_12 = 1._wp / 12._wp
zdt = berg_dt
z1_dt = 1._wp / zdt
!
! we're either going to ignore berg fresh water melt flux and associated heat
! or we pass it into the ocean, so at this point we set them both to zero,
! accumulate the contributions to them from each iceberg in the while loop following
! and then pass them (or not) to the ocean
!
berg_grid%floating_melt(:,:) = 0._wp
! calving_hflx re-used here as temporary workspace for the heat flux associated with melting
berg_grid%calving_hflx(:,:) = 0._wp
!
this => first_berg
DO WHILE( ASSOCIATED(this) )
!
pt => this%current_point
nknberg = this%number(1)
CALL icb_utl_interp( pt%xi, pt%yj, & ! position
& pssu=pt%ssu, pua=pt%ua, & ! oce/atm velocities
& pssv=pt%ssv, pva=pt%va, & ! oce/atm velocities
& psst=pt%sst, pcn=pt%cn, &
& psss=pt%sss )
IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN
CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2, &
& pui=pt%ui, pssh_i=pt%ssh_x, &
& pvi=pt%vi, pssh_j=pt%ssh_y, &
& phi=pt%hi, &
& plat=pt%lat, plon=pt%lon )
END IF
!
zSST = pt%sst
zSSS = pt%sss
CALL eos_fzp(zSSS,zfzpt) ! freezing point
zIC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ???
zM = pt%mass
zT = pt%thickness ! total thickness
zD = rho_berg_1_oce * zT ! draught (keel depth)
zW = pt%width
zL = pt%length
zxi = pt%xi ! position in (i,j) referential
zyj = pt%yj
ii = INT( zxi + 0.5 ) ! T-cell of the berg
ii = mi1( ii + (nn_hls-1) )
ij = INT( zyj + 0.5 )
ij = mj1( ij + (nn_hls-1) )
zVol = zT * zW * zL
! Environment
! default sst, ssu and ssv
! ln_M2016: use temp, u and v profile
IF ( ln_M2016 ) THEN
! load t, u, v and e3 profile at icb position
CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t )
!compute bottom level
CALL icb_utl_getkb( pt%kb, ze3t, zD )
ikb = MIN(pt%kb,mbkt(ii,ij)) ! limit pt%kb by mbkt
! => bottom temperature used to fill ztoce(mbkt:jpk)
ztb = ztoce(ikb) ! basal temperature
zub = zuoce(ikb)
zvb = zvoce(ikb)
ELSE
ztb = pt%sst
zub = pt%ssu
zvb = pt%ssv
END IF
zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 ) ! relative basal velocity
zdva = SQRT( (pt%ua -pt%ssu)**2 + (pt%va -pt%ssv)**2 ) ! relative wind
zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9)
!
! Melt rates in m/s (i.e. division by rday)
! Buoyant convection at sides (eqn M.A10)
IF ( ln_M2016 ) THEN
! averaging along all the iceberg draft
zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday
CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb )
ELSE
zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday
END IF
!
! Basal turbulent melting (eqn M.A7 )
IF ( zSST > zfzpt ) THEN ! Calculate basal melting only if SST above freezing point
zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday
ELSE
zMb = 0._wp ! No basal melting if SST below freezing point
ENDIF
!
! Wave erosion (eqn M.A8 )
zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday
IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass
zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m)
znVol = zTn * zW * zL ! new volume (m^3)
zMnew1 = ( znVol / zVol ) * zM ! new mass (kg)
zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg)
!
zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m)
zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m)
znVol = zTn * zWn * zLn ! new volume (m^3)
zMnew2 = ( znVol / zVol ) * zM ! new mass (kg)
zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg)
!
zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m)
zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m)
znVol = zTn * zWn * zLn ! new volume (m^3)
zMnew = ( znVol / zVol ) * zM ! new mass (kg)
zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg)
zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg)
!
ELSE ! Update dimensions of berg
zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m)
zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m)
zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m)
! Update volume and mass of berg
znVol = zTn*zWn*zLn ! (m^3)
zMnew = (znVol/zVol)*zM ! (kg)
zdM = zM - zMnew ! (kg)
zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg)
zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg)
zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg)
ENDIF
IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits
!
zMbits = pt%mass_of_bits ! mass of bergy bits (kg)
zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg)
znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg)
zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters
zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits)
zMbb = MAX( 0.58_wp*(zdvob**0.8_wp)*(zSST+2._wp) / &
& ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits)
zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s
zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg)
znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt
IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then
zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits
znMbits = 0._wp
ENDIF
ELSE ! No bergy bits
zAbits = 0._wp
zdMbitsE = 0._wp
zdMbitsM = 0._wp
znMbits = pt%mass_of_bits ! retain previous value incase non-zero
ENDIF
! use tmask rather than tmask_i when dealing with icebergs
IF( tmask(ii,ij,1) /= 0._wp ) THEN ! Add melting to the grid and field diagnostics
z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
z1_dt_e1e2 = z1_dt * z1_e1e2
!
! iceberg melt
!! the use of DDPDD function for the cumulative sum is needed for reproducibility
zmelt = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt ! kg/s
CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, dp ), cicb_melt(ii,ij) )
!
! iceberg heat flux
!! the use of DDPDD function for the cumulative sum is needed for reproducibility
!! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the
!! heat density of the icebergs is zero and the heat content flux to the ocean from iceberg
!! melting is always zero. Leaving the term in the code until such a time as this is fixed. DS.
zheat_hcflux = zmelt * pt%heat_density ! heat content flux : kg/s x J/kg = J/s
zheat_latent = - zmelt * rLfus ! latent heat flux: kg/s x J/kg = J/s
CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, dp ), cicb_hflx(ii,ij) )
!
! diagnostics
CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling, &
& zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, &
& zdMv, z1_dt_e1e2, z1_e1e2 )
ELSE
WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij
CALL icb_utl_print_berg( this, kt )
WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij)
CALL ctl_stop('icb_thm', 'berg appears to have grounded!')
ENDIF
! Rolling
zDn = rho_berg_1_oce * zTn ! draught (keel depth)
IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN
zT = zTn
zTn = zWn
zWn = zT
ENDIF
! Store the new state of iceberg (with L>W)
pt%mass = zMnew
pt%mass_of_bits = znMbits
pt%thickness = zTn
pt%width = MIN( zWn , zLn )
pt%length = MAX( zWn , zLn )
next=>this%next
!!gm add a test to avoid over melting ?
!!pm I agree, over melting could break conservation (more melt than calving)
IF( zMnew <= 0._wp ) THEN ! Delete the berg if completely melted
CALL icb_utl_delete( first_berg, this )
!
ELSE ! Diagnose mass distribution on grid
z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling
CALL icb_dia_size( ii, ij, zWn, zLn, zAbits, &
& this%mass_scaling, zMnew, znMbits, z1_e1e2 )
ENDIF
!
this=>next
!
END DO
!
berg_grid%floating_melt = REAL(cicb_melt,dp) ! kg/m2/s
berg_grid%calving_hflx = REAL(cicb_hflx,dp)
!
! now use melt and associated heat flux in ocean (or not)
!
IF(.NOT. ln_passive_mode ) THEN
emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:)
qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:)
ENDIF
!
END SUBROUTINE icb_thm
!!======================================================================
END MODULE icbthm