Newer
Older
MODULE isfcavmlt
!!======================================================================
!! *** MODULE isfcavmlt ***
!! ice shelf module : update surface ocean boundary condition under ice
!! shelves
!!======================================================================
!! History : 4.0 ! 2019-09 (P. Mathiot) Original code
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! isfcav_mlt : compute or read ice shelf fwf/heat fluxes from isf
!! to oce
!!----------------------------------------------------------------------
USE isf_oce ! ice shelf
USE isftbl , ONLY: isf_tbl ! ice shelf depth average
USE isfutils,ONLY: debug ! debug subroutine
USE dom_oce ! ocean space and time domain
USE phycst , ONLY: rcp, rho0, rho0_rcp ! physical constants
USE eosbn2 , ONLY: eos_fzp ! equation of state
USE in_out_manager ! I/O manager
USE iom , ONLY: iom_put ! I/O library
USE fldread , ONLY: fld_read, FLD, FLD_N !
USE lib_fortran, ONLY: glob_sum !
USE lib_mpp , ONLY: ctl_stop !
IMPLICIT NONE
PRIVATE
PUBLIC isfcav_mlt
!! * Substitutions
sparonuz
committed
# include "single_precision_substitute.h90"
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
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
! -------------------------------------------------------------------------------------------------------
! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
! -------------------------------------------------------------------------------------------------------
SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl, &
& pqhc, pqoce, pqfwf )
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt ***
!!
!! ** Purpose : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
!!
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat and fwf fluxes
!!-------------------------- IN -------------------------------------
INTEGER, INTENT(in) :: kt
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pgt , pgs ! gamma t and gamma s
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! top boundary layer tracer
!!---------------------------------------------------------------------
!
! compute latent heat and melt (2d)
SELECT CASE ( cn_isfcav_mlt )
CASE ( 'spe' ) ! ice shelf melt specified (read input file, and heat fluxes derived from
CALL isfcav_mlt_spe( kt, pstbl, &
& pqhc, pqoce, pqfwf )
CASE ( '2eq' ) ! ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
CALL isfcav_mlt_2eq( pgt, pttbl, pstbl, &
& pqhc , pqoce, pqfwf )
CASE ( '3eq' ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
CALL isfcav_mlt_3eq( pgt, pgs , pttbl, pstbl, &
& pqhc, pqoce, pqfwf )
CASE ( 'oasis' ) ! fwf pass trough oasis
CALL isfcav_mlt_oasis( kt, pstbl, &
& pqhc, pqoce, pqfwf )
CASE DEFAULT
CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
END SELECT
!
IF (ln_isfdebug) THEN
IF(lwp) WRITE(numout,*) ''
CALL debug( 'isfcav_mlt qhc :', pqhc (:,:) )
CALL debug( 'isfcav_mlt qoce :', pqoce(:,:) )
CALL debug( 'isfcav_mlt qfwf :', pqfwf(:,:) )
IF(lwp) WRITE(numout,*) ''
END IF
!
END SUBROUTINE isfcav_mlt
! -------------------------------------------------------------------------------------------------------
! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
! -------------------------------------------------------------------------------------------------------
SUBROUTINE isfcav_mlt_spe(kt, pstbl, & ! <<== in
& pqhc , pqoce, pqfwf ) ! ==>> out
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt_spe ***
!!
!! ** Purpose : - read ice shelf melt from forcing file
!! - compute ocea-ice heat flux (assuming it is equal to latent heat)
!! - compute heat content flux
!!---------------------------------------------------------------------
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat content, latent heat and fwf fluxes
!!-------------------------- IN -------------------------------------
INTEGER , INTENT(in ) :: kt ! current time step
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pstbl ! salinity in tbl
!!--------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! tbl freezing temperature
!!--------------------------------------------------------------------
!
! Compute freezing temperature
sparonuz
committed
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), CASTDP(risfdep(:,:)) )
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
!
! read input file of fwf (from isf to oce; ie melt)
CALL fld_read ( kt, 1, sf_isfcav_fwf )
!
! define fwf and qoce
! ocean heat flux is assume to be equal to the latent heat
pqfwf(:,:) = sf_isfcav_fwf(1)%fnow(:,:,1) ! fwf ( > 0 from isf to oce)
pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean heat flux ( > 0 from isf to oce)
pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
!
! output freezing point at the interface
CALL iom_put('isftfrz_cav', ztfrz(:,:) * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_spe
SUBROUTINE isfcav_mlt_2eq(pgt , pttbl, pstbl, & ! <<== in
& pqhc, pqoce, pqfwf ) ! ==>> out
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt_2eq ***
!!
!! ** Purpose : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
!!
!! ** Method : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux.
!! From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being :
!! qfwf = Gammat * rho0 * Cp * ( Tw - Tfrz ) / Lf
!! qhoce = qlat
!! qhc = qfwf * Cp * Tfrz
!!
!! ** Reference : Hunter, J. R.: Specification for test models of ice shelf cavities,
!! Tech. Rep. June, Antarctic Climate & Ecosystems Cooperative Research Centre, available at:
!! http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
!!
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! hean content, ocean-ice heat and fwf fluxes
!!-------------------------- IN -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pgt ! temperature exchange coeficient
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! temperature and salinity in top boundary layer
!!--------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! freezing temperature
REAL(wp), DIMENSION(jpi,jpj) :: zthd ! thermal driving
!!--------------------------------------------------------------------
!
! Calculate freezing temperature
sparonuz
committed
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), CASTDP(risfdep(:,:)) )
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
!
! thermal driving
zthd (:,:) = ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:)
!
! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
pqfwf(:,:) = pgt(:,:) * rho0_rcp * zthd(:,:) / rLfusisf ! fresh water flux ( > 0 from isf to oce)
pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocea-ice flux ( > 0 from isf to oce)
pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
!
! output thermal driving and freezinpoint at the ice shelf interface
CALL iom_put('isfthermald_cav', zthd )
CALL iom_put('isftfrz_cav' , ztfrz(:,:) * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_2eq
SUBROUTINE isfcav_mlt_3eq(pgt, pgs , pttbl, pstbl, & ! <<== in
& pqhc, pqoce, pqfwf ) ! ==>> out
!!----------------------------------------------------------------------
!!
!! *** ROUTINE isfcav_mlt_3eq ***
!!
!! ** Purpose : Compute ice shelf fwf/heqt fluxes using the 3 equation formulation
!!
!! ** Method : The melt rate is determined considering the heat balance, the salt balance
!! at the phase change interface and a linearisation of the equation of state.
!!
!! ** Reference : - Holland, D. M. and Jenkins, A.,
!! Modeling Thermodynamic Ice-Ocean Interactions at the Base of an Ice Shelf,
!! J. Phys. Oceanogr., 29, 1999.
!! - Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone,
!! R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland,
!! P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.:
!! Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects:
!! MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1),
!! Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016.
!!
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! latent heat and fwf fluxes
!!-------------------------- IN -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pgt , pgs ! heat/salt exchange coeficient
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pttbl, pstbl ! mean temperature and salinity in top boundary layer
!!--------------------------------------------------------------------
REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 ! dummy local scalar for quadratic equation resolution
REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac ! dummy local scalar for quadratic equation resolution
REAL(wp) :: zeps = 1.e-20
REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! freezing point
REAL(wp), DIMENSION(jpi,jpj) :: zqcon ! conductive flux through the ice shelf
REAL(wp), DIMENSION(jpi,jpj) :: zthd ! thermal driving
!
INTEGER :: ji, jj ! dummy loop indices
!!--------------------------------------------------------------------
!
! compute upward heat flux zhtflx and upward water flux zwflx
! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
! compute coeficient to solve the 2nd order equation
zeps1 = rho0_rcp * pgt(ji,jj)
zeps2 = rLfusisf * rho0 * pgs(ji,jj)
zeps3 = rhoisf * rcpisf * rkappa / MAX(risfdep(ji,jj),zeps)
zeps4 = risf_lamb2 + risf_lamb3 * risfdep(ji,jj)
zeps6 = zeps4 - pttbl(ji,jj)
zeps7 = zeps4 - rtsurf
!
! solve the 2nd order equation to find zsfrz
zaqe = risf_lamb1 * (zeps1 + zeps3)
zaqer = 0.5_wp / MIN(zaqe,-zeps)
zbqe = zeps1 * zeps6 + zeps3 * zeps7 - zeps2
zcqe = zeps2 * pstbl(ji,jj)
zdis = zbqe * zbqe - 4.0_wp * zaqe * zcqe
!
! Presumably zdis can never be negative because gammas is very small compared to gammat
zsfrz=(-zbqe - SQRT(zdis)) * zaqer
IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe + SQRT(zdis)) * zaqer ! check this if this if is needed
!
! compute t freeze (eq. 25)
ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
!
! thermal driving
zthd(ji,jj) = ( pttbl(ji,jj) - ztfrz(ji,jj) )
!
! compute the upward water and heat flux (eq. 24 and eq. 26)
pqfwf(ji,jj) = - rho0 * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux ( > 0 from isf to oce)
pqoce(ji,jj) = - rho0_rcp * pgt(ji,jj) * zthd (ji,jj) ! ocean-ice heat flux ( > 0 from isf to oce)
pqhc (ji,jj) = rcp * pqfwf(ji,jj) * ztfrz(ji,jj) ! heat content flux ( > 0 from isf to oce)
!
zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )
!
END_2D
!
! output conductive heat flux through the ice
CALL iom_put('qconisf', zqcon(:,:) * mskisf_cav(:,:) )
!
! output thermal driving and freezing point at the interface
CALL iom_put('isfthermald_cav', zthd (:,:) * mskisf_cav(:,:) )
CALL iom_put('isftfrz_cav' , ztfrz(:,:) * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_3eq
SUBROUTINE isfcav_mlt_oasis(kt, pstbl, & ! <<== in
& pqhc , pqoce, pqfwf ) ! ==>> out
!!----------------------------------------------------------------------
!! *** ROUTINE isfcav_mlt_oasis ***
!!
!! ** Purpose : scale the fwf read from input file by the total amount received by the sbccpl interface
!!
!! ** Purpose : - read ice shelf melt from forcing file => pattern
!! - total amount of fwf is given by sbccpl (fwfisf_cpl)
!! - scale fwf and compute heat fluxes
!!
!!---------------------------------------------------------------------
!!-------------------------- OUT -------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pqhc, pqoce, pqfwf ! heat content, latent heat and fwf fluxes
!!-------------------------- IN -------------------------------------
INTEGER , INTENT(in ) :: kt ! current time step
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pstbl ! salinity in tbl
!!--------------------------------------------------------------------
REAL(wp) :: zfwf_fld, zfwf_oasis ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
REAL(wp), DIMENSION(jpi,jpj) :: ztfrz ! tbl freezing temperature
REAL(wp), DIMENSION(jpi,jpj) :: zfwf ! 2d fwf map after scaling
!!--------------------------------------------------------------------
!
! Calculate freezing temperature
sparonuz
committed
CALL eos_fzp( pstbl(:,:), ztfrz(:,:), CASTDP(risfdep(:,:)) )
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
!
! read input file of fwf from isf to oce
CALL fld_read ( kt, 1, sf_isfcav_fwf )
!
! ice shelf 2d map
zfwf(:,:) = sf_isfcav_fwf(1)%fnow(:,:,1)
!
! compute glob sum from input file
! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
zfwf_fld = glob_sum('isfcav_mlt', e1e2t(:,:) * zfwf(:,:))
!
! compute glob sum from atm->oce ice shelf fwf
! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(:,:) * fwfisf_oasis(:,:))
!
! scale fwf
zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
!
! define fwf and qoce
! ocean heat flux is assume to be equal to the latent heat
pqfwf(:,:) = zfwf(:,:) ! fwf ( > 0 from isf to oce)
pqoce(:,:) = - pqfwf(:,:) * rLfusisf ! ocean heat flux ( > 0 from isf to oce)
pqhc (:,:) = pqfwf(:,:) * ztfrz(:,:) * rcp ! heat content flux ( > 0 from isf to oce)
!
CALL iom_put('isftfrz_cav', ztfrz * mskisf_cav(:,:) )
!
END SUBROUTINE isfcav_mlt_oasis
END MODULE isfcavmlt