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
MODULE icevar
!!======================================================================
!! *** MODULE icevar ***
!! sea-ice: series of functions to transform or compute ice variables
!!======================================================================
!! History : - ! 2006-01 (M. Vancoppenolle) Original code
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!----------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!!
!! There are three sets of variables
!! VGLO : global variables of the model
!! - v_i (jpi,jpj,jpl)
!! - v_s (jpi,jpj,jpl)
!! - a_i (jpi,jpj,jpl)
!! - t_s (jpi,jpj,jpl)
!! - e_i (jpi,jpj,nlay_i,jpl)
!! - e_s (jpi,jpj,nlay_s,jpl)
!! - sv_i(jpi,jpj,jpl)
!! - oa_i(jpi,jpj,jpl)
!! VEQV : equivalent variables sometimes used in the model
!! - h_i(jpi,jpj,jpl)
!! - h_s(jpi,jpj,jpl)
!! - t_i(jpi,jpj,nlay_i,jpl)
!! ...
!! VAGG : aggregate variables, averaged/summed over all
!! thickness categories
!! - vt_i(jpi,jpj)
!! - vt_s(jpi,jpj)
!! - at_i(jpi,jpj)
!! - st_i(jpi,jpj)
!! - et_s(jpi,jpj) total snow heat content
!! - et_i(jpi,jpj) total ice thermal content
!! - sm_i(jpi,jpj) mean ice salinity
!! - tm_i(jpi,jpj) mean ice temperature
!! - tm_s(jpi,jpj) mean snw temperature
!!----------------------------------------------------------------------
!! ice_var_agg : integrate variables over layers and categories
!! ice_var_glo2eqv : transform from VGLO to VEQV
!! ice_var_eqv2glo : transform from VEQV to VGLO
!! ice_var_salprof : salinity profile in the ice
!! ice_var_salprof1d : salinity profile in the ice 1D
!! ice_var_zapsmall : remove very small area and volume
!! ice_var_zapneg : remove negative ice fields
!! ice_var_roundoff : remove negative values arising from roundoff erros
!! ice_var_brine : brine volume
!! ice_var_enthalpy : compute ice and snow enthalpies from temperature
!! ice_var_sshdyn : compute equivalent ssh in lead
!! ice_var_itd : convert N-cat to M-cat
!! ice_var_snwfra : fraction of ice covered by snow
!! ice_var_snwblow : distribute snow fall between ice and ocean
!!----------------------------------------------------------------------
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants (ocean directory)
USE sbc_oce , ONLY : sss_m, sst_m, ln_ice_embd, nn_fsbc
USE ice ! sea-ice: variables
USE ice1D ! sea-ice: thermodynamics variables
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
IMPLICIT NONE
PRIVATE
PUBLIC ice_var_agg
PUBLIC ice_var_glo2eqv
PUBLIC ice_var_salprof
PUBLIC ice_var_salprof1d
PUBLIC ice_var_zapsmall
PUBLIC ice_var_zapneg
PUBLIC ice_var_roundoff
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
PUBLIC ice_var_sshdyn
PUBLIC ice_var_itd
PUBLIC ice_var_snwfra
PUBLIC ice_var_snwblow
INTERFACE ice_var_itd
MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
END INTERFACE
!! * Substitutions
# include "do_loop_substitute.h90"
INTERFACE ice_var_snwfra
MODULE PROCEDURE ice_var_snwfra_1d, ice_var_snwfra_2d, ice_var_snwfra_3d
END INTERFACE
INTERFACE ice_var_snwblow
MODULE PROCEDURE ice_var_snwblow_1d, ice_var_snwblow_2d
END INTERFACE
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icevar.F90 15385 2021-10-15 13:52:48Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_var_agg( kn )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_agg ***
!!
!! ** Purpose : aggregates ice-thickness-category variables to
!! all-ice variables, i.e. it turns VGLO into VAGG
!!-------------------------------------------------------------------
INTEGER, INTENT( in ) :: kn ! =1 state variables only
! ! >1 state variables + others
!
INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: z1_vt_i, z1_vt_s, z1_at_i
!!-------------------------------------------------------------------
!
! full arrays: vt_i, vt_s, at_i, vt_ip, vt_il, at_ip
! reduced arrays: the rest
! --- integrated values
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
vt_i(ji,jj) = SUM( v_i (ji,jj,:) )
vt_s(ji,jj) = SUM( v_s (ji,jj,:) )
at_i(ji,jj) = SUM( a_i (ji,jj,:) )
!
at_ip(ji,jj) = SUM( a_ip(ji,jj,:) ) ! melt ponds
vt_ip(ji,jj) = SUM( v_ip(ji,jj,:) )
vt_il(ji,jj) = SUM( v_il(ji,jj,:) )
!
ato_i(ji,jj) = 1._wp - at_i(ji,jj) ! open water fraction
END_2D
DO_2D( 0, 0, 0, 0 )
et_s(ji,jj) = SUM( SUM( e_s (ji,jj,:,:), dim=2 ) )
et_i(ji,jj) = SUM( SUM( e_i (ji,jj,:,:), dim=2 ) )
!
!!GS: tm_su always needed by ABL over sea-ice
IF( at_i(ji,jj) <= epsi20 ) THEN
tm_su (ji,jj) = rt0
ELSE
tm_su (ji,jj) = SUM( t_su(ji,jj,:) * a_i(ji,jj,:) ) / at_i(ji,jj)
ENDIF
END_2D
IF( nn_icesal == 4 ) THEN
DO_2D( 0, 0, 0, 0 )
st_i(ji,jj) = SUM( SUM( szv_i (ji,jj,:,:), dim=2 ) )
END_2D
ELSE
DO_2D( 0, 0, 0, 0 )
st_i(ji,jj) = SUM( sv_i(ji,jj,:) )
END_2D
ENDIF
! The following fields are calculated for diagnostics and outputs only
! ==> Do not use them for other purposes
IF( kn > 1 ) THEN
!
DO_2D( 0, 0, 0, 0 )
IF( at_i(ji,jj) > epsi20 ) THEN ; z1_at_i = 1._wp / at_i(ji,jj)
ELSE ; z1_at_i = 0._wp
ENDIF
IF( vt_i(ji,jj) > epsi20 ) THEN ; z1_vt_i = 1._wp / vt_i(ji,jj)
ELSE ; z1_vt_i = 0._wp
ENDIF
! mean ice/snow thickness
hm_i(ji,jj) = vt_i(ji,jj) * z1_at_i
hm_s(ji,jj) = vt_s(ji,jj) * z1_at_i
!
! mean temperature (K), salinity and age
tm_si(ji,jj) = SUM( t_si(ji,jj,:) * a_i(ji,jj,:) ) * z1_at_i
om_i (ji,jj) = SUM( oa_i(ji,jj,:) ) * z1_at_i
sm_i (ji,jj) = st_i(ji,jj) * z1_vt_i
END_2D
!
tm_i(:,:) = 0._wp
tm_s(:,:) = 0._wp
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_i )
IF( vt_i(ji,jj) > epsi20 ) THEN
tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * t_i (ji,jj,jk,jl) * v_i(ji,jj,jl) / vt_i(ji,jj)
ELSE
tm_i(ji,jj) = rt0
ENDIF
END_3D
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
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_s )
IF( vt_s(ji,jj) > epsi20 ) THEN
tm_s(ji,jj) = tm_s(ji,jj) + r1_nlay_s * t_s (ji,jj,jk,jl) * v_s(ji,jj,jl) / vt_s(ji,jj)
ELSE
tm_s(ji,jj) = rt0
ENDIF
END_3D
END DO
!
!!$ DO_2D( 0, 0, 0, 0 )
!!$ IF( at_i(ji,jj) > epsi20 ) THEN ; z1_at_i = 1._wp / at_i(ji,jj)
!!$ ELSE ; z1_at_i = 0._wp
!!$ ENDIF
!!$ IF( vt_i(ji,jj) > epsi20 ) THEN ; z1_vt_i = 1._wp / vt_i(ji,jj)
!!$ ELSE ; z1_vt_i = 0._wp
!!$ ENDIF
!!$ IF( vt_s(ji,jj) > epsi20 ) THEN ; z1_vt_s = 1._wp / vt_s(ji,jj)
!!$ ELSE ; z1_vt_s = 0._wp
!!$ ENDIF
!!$
!!$ ! mean ice/snow thickness
!!$ hm_i(ji,jj) = vt_i(ji,jj) * z1_at_i
!!$ hm_s(ji,jj) = vt_s(ji,jj) * z1_at_i
!!$ !
!!$ ! mean temperature (K), salinity and age
!!$ tm_si(ji,jj) = SUM( t_si(ji,jj,:) * a_i(ji,jj,:) ) * z1_at_i
!!$ om_i (ji,jj) = SUM( oa_i(ji,jj,:) ) * z1_at_i
!!$ sm_i (ji,jj) = st_i(ji,jj) * z1_vt_i
!!$ !
!!$ tm_i(ji,jj) = 0._wp
!!$ tm_s(ji,jj) = 0._wp
!!$ DO jl = 1, jpl
!!$ DO jk = 1, nlay_i
!!$ tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * t_i (ji,jj,jk,jl) * v_i(ji,jj,jl) * z1_vt_i
!!$ END DO
!!$ DO jk = 1, nlay_s
!!$ tm_s(ji,jj) = tm_s(ji,jj) + r1_nlay_s * t_s (ji,jj,jk,jl) * v_s(ji,jj,jl) * z1_vt_s
!!$ END DO
!!$ END DO
!!$ !
!!$ END_2D
! put rt0 where there is no ice
WHERE( at_i(A2D(0)) <= epsi20 )
tm_si(:,:) = rt0
tm_i (:,:) = rt0
tm_s (:,:) = rt0
END WHERE
!
! mean melt pond depth
WHERE( at_ip(A2D(0)) > epsi20 )
hm_ip(:,:) = vt_ip(A2D(0)) / at_ip(A2D(0))
hm_il(:,:) = vt_il(A2D(0)) / at_ip(A2D(0))
ELSEWHERE
hm_ip(:,:) = 0._wp
hm_il(:,:) = 0._wp
END WHERE
!
ENDIF
!
END SUBROUTINE ice_var_agg
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_glo2eqv ***
!!
!! ** Purpose : computes equivalent variables as function of
!! global variables, i.e. it turns VGLO into VEQV
!!-------------------------------------------------------------------
INTEGER, INTENT( in ) :: kn ! =1 everything including ponds (necessary for init)
! ! =2 excluding ponds if ln_pnd=F
REAL(wp) :: ze_i, zs_i ! local scalars
REAL(wp) :: zhmax, z1_hmax ! - -
REAL(wp) :: zlay_i, zlay_s ! - -
REAL(wp), PARAMETER :: zhl_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation
REAL(wp), PARAMETER :: zhl_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation
REAL(wp) :: z1_hl, z1_a_i, z1_a_ip
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: za_s_fra
!!-------------------------------------------------------------------
!!gm Question 2: It is possible to define existence of sea-ice in a common way between
!! ice area and ice volume ?
!! the idea is to be able to define one for all at the begining of this routine
!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
!---------------------------------------------------------------
! Ice thickness, snow thickness, ice salinity, ice age and ponds
!---------------------------------------------------------------
!
! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
! clem: if a>1 then do something
zhmax = hi_max(jpl)
z1_hmax = 1._wp / hi_max(jpl)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( v_i(ji,jj,jpl) > ( zhmax*a_i(ji,jj,jpl) ) ) a_i(ji,jj,jpl) = MIN( 1._wp, v_i(ji,jj,jpl) * z1_hmax )
END_2D
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! !--- inverse of the ice area
IF( a_i(ji,jj,jl) > epsi20 ) THEN ; z1_a_i = 1._wp / a_i(ji,jj,jl)
ELSE ; z1_a_i = 0._wp
ENDIF
! !--- ice thickness
h_i(ji,jj,jl) = v_i (ji,jj,jl) * z1_a_i
! !--- snow thickness
h_s(ji,jj,jl) = v_s (ji,jj,jl) * z1_a_i
! !--- ice age
o_i(ji,jj,jl) = oa_i(ji,jj,jl) * z1_a_i
!
END_2D
ENDDO
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
IF( kn == 1 .OR. ln_pnd ) THEN
ALLOCATE( za_s_fra(A2D(0),jpl) )
!
z1_hl = 1._wp / ( zhl_max - zhl_min )
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( a_ip(ji,jj,jl) > epsi20 ) THEN ; z1_a_ip = 1._wp / a_ip(ji,jj,jl)
ELSE ; z1_a_ip = 0._wp
ENDIF
! !--- pond and lid thickness
h_ip(ji,jj,jl) = v_ip(ji,jj,jl) * z1_a_ip
h_il(ji,jj,jl) = v_il(ji,jj,jl) * z1_a_ip
END_2D
! !--- melt pond effective area (used for albedo)
DO_2D( 0, 0, 0, 0 )
IF( a_i(ji,jj,jl) > epsi20 ) THEN ; a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl)
ELSE ; a_ip_frac(ji,jj,jl) = 0._wp
ENDIF
IF ( h_il(ji,jj,jl) <= zhl_min ) THEN ; a_ip_eff(ji,jj,jl) = a_ip_frac(ji,jj,jl) ! lid is very thin. Expose all the pond
ELSEIF( h_il(ji,jj,jl) >= zhl_max ) THEN ; a_ip_eff(ji,jj,jl) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow
ELSE ; a_ip_eff(ji,jj,jl) = a_ip_frac(ji,jj,jl) * & ! lid is in between. Expose part of the pond
& ( zhl_max - h_il(ji,jj,jl) ) * z1_hl
ENDIF
!
END_2D
ENDDO
!
CALL ice_var_snwfra( h_s(A2D(0),:), za_s_fra(:,:,:) ) ! calculate ice fraction covered by snow
a_ip_eff(:,:,:) = MIN( a_ip_eff(:,:,:), 1._wp - za_s_fra(:,:,:) ) ! make sure (a_ip_eff + a_s_fra) <= 1
!
DEALLOCATE( za_s_fra )
ENDIF
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
!
!-------------------
! Ice salinity (with a min value rn_simin and a max value rn_sinew*sss)
!-------------------
IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN
!
CALL ice_var_salprof ! salinity profile
!
ELSEIF( nn_icesal == 2 ) THEN
!
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF( v_i(ji,jj,jl) > epsi20 ) THEN
!!clem test s_i(ji,jj,jl) = MAX( rn_simin , MIN( rn_sinew * sss_m(ji,jj), sv_i(ji,jj,jl) / v_i(ji,jj,jl) ) )
s_i(ji,jj,jl) = sv_i(ji,jj,jl) / v_i(ji,jj,jl)
ELSE
s_i(ji,jj,jl) = rn_simin
ENDIF
END_2D
ENDDO
CALL ice_var_salprof ! salinity profile
!
ELSEIF( nn_icesal == 4 ) THEN
!
s_i(:,:,:) = 0._wp
zlay_i = REAL( nlay_i , wp ) ! number of layers
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
zs_i = szv_i(ji,jj,jk,jl) / ( v_i(ji,jj,jl) * r1_nlay_i )
!!clem test sz_i(ji,jj,jk,jl) = MAX( rn_simin , MIN( rn_sinew * sss_m(ji,jj), zs_i ) )
sz_i(ji,jj,jk,jl) = zs_i
ELSE !--- no ice
sz_i(ji,jj,jk,jl) = rn_simin
ENDIF
s_i(ji,jj,jl) = s_i(ji,jj,jl) + sz_i(ji,jj,jk,jl) * r1_nlay_i
END_3D
END DO
!
ENDIF
!-------------------
! Ice temperature [K] (with a minimum value (rt0 - 100.))
!-------------------
zlay_i = REAL( nlay_i , wp ) ! number of layers
DO jk = 1, nlay_i
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
!
ze_i = e_i (ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
! Conversion q(S,T) -> T (second order equation)
zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
!
ELSE !--- no ice
t_i(ji,jj,jk,jl) = rt0
ENDIF
END_2D
END DO
END DO
!--------------------
! Snow temperature [K] (with a minimum value (rt0 - 100.))
!--------------------
zlay_s = REAL( nlay_s , wp )
DO jl = 1, jpl
DO jk = 1, nlay_s
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( v_s(ji,jj,jl) > epsi20 ) THEN !--- icy area
t_s(ji,jj,jk,jl) = rt0 + MAX( -100._wp , &
& MIN( r1_rcpi*( -r1_rhos*( e_s(ji,jj,jk,jl) / v_s(ji,jj,jl) * zlay_s ) + rLfus ) , 0._wp ) )
ELSE !--- no ice
t_s(ji,jj,jk,jl) = rt0
ENDIF
END_2D
END DO
vt_i (:,:) = SUM( v_i, dim=3 )
vt_s (:,:) = SUM( v_s, dim=3 )
at_i (:,:) = SUM( a_i, dim=3 )
!!$ SUBROUTINE ice_var_eqv2glo
!!$ !!-------------------------------------------------------------------
!!$ !! *** ROUTINE ice_var_eqv2glo ***
!!$ !!
!!$ !! ** Purpose : computes global variables as function of
!!$ !! equivalent variables, i.e. it turns VEQV into VGLO
!!$ !!-------------------------------------------------------------------
!!$ INTEGER :: ji, jj, jl ! dummy loop indices
!!$ !!-------------------------------------------------------------------
!!$ !
!!$ DO jl = 1, jpl
!!$ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!!$ v_i (ji,jj,jl) = h_i (ji,jj,jl) * a_i (ji,jj,jl)
!!$ v_s (ji,jj,jl) = h_s (ji,jj,jl) * a_i (ji,jj,jl)
!!$ sv_i(ji,jj,jl) = s_i (ji,jj,jl) * v_i (ji,jj,jl)
!!$ v_ip(ji,jj,jl) = h_ip(ji,jj,jl) * a_ip(ji,jj,jl)
!!$ v_il(ji,jj,jl) = h_il(ji,jj,jl) * a_ip(ji,jj,jl)
!!$ END_2D
!!$ ENDDO
!!$ !
!!$ END SUBROUTINE ice_var_eqv2glo
SUBROUTINE ice_var_salprof
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_salprof ***
!!
!! ** Purpose : computes salinity profile in function of bulk salinity
!!
!! ** Method : If bulk salinity greater than zsi1,
!! the profile is assumed to be constant (S_inf)
!! If bulk salinity lower than zsi0,
!! the profile is linear with 0 at the surface (S_zero)
!! If it is between zsi0 and zsi1, it is a
!! alpha-weighted linear combination of s_inf and s_zero
!!
!! ** References : Vancoppenolle et al., 2007
!!-------------------------------------------------------------------
INTEGER :: ji, jj, jk, jl ! dummy loop index
REAL(wp) :: z1_dS, ztmp1, ztmp2, zalpha
REAL(wp), PARAMETER :: zsi0 = 3.5_wp
REAL(wp), PARAMETER :: zsi1 = 4.5_wp
!!-------------------------------------------------------------------
!!gm Question: Remove the option 3 ? How many years since it last use ?
SELECT CASE ( nn_icesal )
!
! !---------------------------------------!
CASE( 1 ) ! constant salinity in time and space !
! !---------------------------------------!
sz_i(:,:,:,:) = rn_icesal
s_i (:,:,:) = rn_icesal
!
! !---------------------------------------------!
CASE( 2 , 4 ) ! time varying salinity with linear profile !
! !---------------------------------------------!
z1_dS = 1._wp / ( zsi1 - zsi0 )
!
DO jl = 1, jpl
DO jk = 1, nlay_i
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
!
IF ( s_i(ji,jj,jl) >= zsi1 ) THEN ; zalpha = 0._wp
ELSEIF( s_i(ji,jj,jl) <= zsi0 ) THEN ; zalpha = 1._wp
ELSE ; zalpha = ( zsi1 - s_i(ji,jj,jl) ) * z1_dS
IF( s_i(ji,jj,jl) >= ( 0.5_wp * rn_sinew * sss_m(ji,jj) ) ) zalpha = 0._wp ! force constant profile when SSS too low (Baltic Sea)
IF( s_i(ji,jj,jl) <= ( REAL( nlay_i , wp ) * rn_simin ) ) zalpha = 0._wp ! force constant profile when ice surface salinity too small
! ! it depends on nlay_i which is not ideal
! linear profile with 0 surface value
sz_i(ji,jj,jk,jl) = zalpha * s_i(ji,jj,jl) * 2._wp * ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i &
& + ( 1._wp - zalpha ) * s_i(ji,jj,jl)
! => mean(sz_i(jk)) = s_i
!
! !-------------------------------------------!
CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
! !-------------------------------------------! (mean = 2.30)
!
s_i(:,:,:) = 2.30_wp
!!gm Remark: if we keep the case 3, then compute an store one for all time-step
!! a array S_prof(1:nlay_i) containing the calculation and just do:
! DO jk = 1, nlay_i
! sz_i(:,:,jk,:) = S_prof(jk)
! END DO
!!gm end
DO jl = 1, jpl
DO jk = 1, nlay_i
ztmp1 = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
ztmp2 = 1.6_wp * ( 1._wp - COS( rpi * ztmp1**(0.407_wp/(0.573_wp+ztmp1)) ) )
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
sz_i(ji,jj,jk,jl) = ztmp2
END_2D
END DO
END DO
END SELECT
!
END SUBROUTINE ice_var_salprof
SUBROUTINE ice_var_salprof1d
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_salprof1d ***
!!
!! ** Purpose : 1d computation of the sea ice salinity profile
!! Works with 1d vectors and is used by thermodynamic modules
!!-------------------------------------------------------------------
INTEGER :: ji, jk ! dummy loop indices
REAL(wp) :: z1_dS, ztmp1, ztmp2, zalpha
REAL(wp), PARAMETER :: zsi0 = 3.5_wp
REAL(wp), PARAMETER :: zsi1 = 4.5_wp
!!-------------------------------------------------------------------
!
SELECT CASE ( nn_icesal )
!
! !---------------------------------------!
CASE( 1 ) ! constant salinity in time and space !
! !---------------------------------------!
sz_i_1d(1:npti,:) = rn_icesal
!
! !---------------------------------------------!
CASE( 2 , 4 ) ! time varying salinity with linear profile !
! !---------------------------------------------!
z1_dS = 1._wp / ( zsi1 - zsi0 )
!
DO jk = 1, nlay_i
DO ji = 1, npti
IF ( s_i_1d(ji) >= zsi1 ) THEN ; zalpha = 0._wp
ELSEIF( s_i_1d(ji) <= zsi0 ) THEN ; zalpha = 1._wp
ELSE ; zalpha = ( zsi1 - s_i_1d(ji) ) * z1_dS
ENDIF
IF( s_i_1d(ji) >= ( 0.5_wp * rn_sinew * sss_1d(ji) ) ) zalpha = 0._wp ! force constant profile when SSS too low (Baltic Sea)
IF( s_i_1d(ji) <= ( REAL( nlay_i , wp ) * rn_simin ) ) zalpha = 0._wp ! force constant profile when ice surface salinity too small
! ! it depends on nlay_i which is not ideal
! linear profile with 0 surface value
sz_i_1d(ji,jk) = zalpha * s_i_1d(ji) * 2._wp * ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i &
& + ( 1._wp - zalpha ) * s_i_1d(ji)
! => mean(sz_i(jk)) = s_i
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
END DO
END DO
!
! !-------------------------------------------!
CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
! !-------------------------------------------! (mean = 2.30)
!
s_i_1d(1:npti) = 2.30_wp
!
!!gm cf remark in ice_var_salprof routine, CASE( 3 )
DO jk = 1, nlay_i
ztmp1 = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
ztmp2 = 1.6_wp * ( 1._wp - COS( rpi * ztmp1**( 0.407_wp / ( 0.573_wp + ztmp1 ) ) ) )
DO ji = 1, npti
sz_i_1d(ji,jk) = ztmp2
END DO
END DO
!
END SELECT
!
END SUBROUTINE ice_var_salprof1d
SUBROUTINE ice_var_zapsmall
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_zapsmall ***
!!
!! ** Purpose : Remove too small sea ice areas and correct fluxes
!!-------------------------------------------------------------------
INTEGER :: ji, jj, jl, jk ! dummy loop indices
!!-------------------------------------------------------------------
!
WHERE( a_i(A2D(0),:) > epsi10 ) ; h_i(A2D(0),:) = v_i(A2D(0),:) / a_i(A2D(0),:)
ELSEWHERE ; h_i(A2D(0),:) = 0._wp
END WHERE
!
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
!-----------------------------------------------------------------
! Zap ice volume, add salt to ocean
!-----------------------------------------------------------------
IF( nn_icesal == 4 ) THEN
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_i )
!
zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) )
!
IF( zsmall < epsi10 ) THEN
! update exchanges with ocean
sfx_res(ji,jj) = sfx_res(ji,jj) + szv_i(ji,jj,jk,jl) * rhoi * r1_Dt_ice
szv_i(ji,jj,jk,jl) = 0._wp
sz_i (ji,jj,jk,jl) = rn_simin
ENDIF
END_3D
ENDDO
ELSE
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
!
zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) )
!
IF( zsmall < epsi10 ) THEN
! update exchanges with ocean
sfx_res(ji,jj) = sfx_res(ji,jj) + sv_i(ji,jj,jl) * rhoi * r1_Dt_ice
sv_i (ji,jj,jl) = 0._wp
ENDIF
END_2D
END DO
ENDIF
!-----------------------------------------------------------------
! Zap ice energy and use ocean heat to melt ice
!-----------------------------------------------------------------
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_i )
!
zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) )
!
IF( zsmall < epsi10 ) THEN
! update exchanges with ocean
hfx_res(ji,jj) = hfx_res(ji,jj) - e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0
e_i(ji,jj,jk,jl) = 0._wp
t_i(ji,jj,jk,jl) = rt0
ENDIF
ENDDO
DO jl = 1, jpl
DO_3D( 0, 0, 0, 0, 1, nlay_s )
zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) )
IF( zsmall < epsi10 ) THEN
! update exchanges with ocean
hfx_res(ji,jj) = hfx_res(ji,jj) - e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0
e_s(ji,jj,jk,jl) = 0._wp
t_s(ji,jj,jk,jl) = rt0
ENDIF
END_3D
ENDDO
!
!-----------------------------------------------------------------
! zap ice and snow volume, add water to ocean
!-----------------------------------------------------------------
DO jl = 1, jpl
DO_2D( 0, 0, 0, 0 )
zsmall = MIN( a_i(ji,jj,jl), v_i(ji,jj,jl), h_i(ji,jj,jl) )
IF( zsmall < epsi10 ) THEN
! update exchanges with ocean
wfx_res(ji,jj) = wfx_res(ji,jj) + v_i (ji,jj,jl) * rhoi * r1_Dt_ice
wfx_res(ji,jj) = wfx_res(ji,jj) + v_s (ji,jj,jl) * rhos * r1_Dt_ice
wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_Dt_ice
!
a_i (ji,jj,jl) = 0._wp
v_i (ji,jj,jl) = 0._wp
v_s (ji,jj,jl) = 0._wp
t_su (ji,jj,jl) = sst_m(ji,jj) + rt0
oa_i (ji,jj,jl) = 0._wp
!
h_i (ji,jj,jl) = 0._wp
h_s (ji,jj,jl) = 0._wp
!
a_ip (ji,jj,jl) = 0._wp
v_ip (ji,jj,jl) = 0._wp
v_il (ji,jj,jl) = 0._wp
h_ip (ji,jj,jl) = 0._wp
h_il (ji,jj,jl) = 0._wp
ENDIF
at_i (A2D(0)) = SUM( a_i (A2D(0),:), dim=3 )
vt_i (A2D(0)) = SUM( v_i (A2D(0),:), dim=3 )
!!clem add?
! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 )
! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 )
! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 )
! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 )
!!clem
! open water = 1 if at_i=0
WHERE( at_i(A2D(0)) == 0._wp ) ato_i(A2D(0)) = 1._wp
SUBROUTINE ice_var_zapneg( ihls, pdt, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i, pszv_i )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_zapneg ***
!!
!! ** Purpose : Remove negative sea ice fields and correct fluxes
!!-------------------------------------------------------------------
INTEGER , INTENT(in ) :: ihls ! loop index
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pszv_i ! ice salt content
!
INTEGER :: ji, jj, jl, jk ! dummy loop indices
REAL(wp) :: z1_dt
REAL(wp), DIMENSION(jpi,jpj) :: zwfx_res, zhfx_res, zsfx_res ! needed since loop is not (0,0,0,0)
!!-------------------------------------------------------------------
!
DO_2D( ihls, ihls, ihls, ihls )
zwfx_res(ji,jj) = 0._wp
zhfx_res(ji,jj) = 0._wp
zsfx_res(ji,jj) = 0._wp
END_2D
! make sure a_i=0 where v_i<=0
WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
!--------------------------------------
! zap ice salt and send it to the ocean
!--------------------------------------
IF( nn_icesal == 4 ) THEN
DO jl = 1, jpl
DO_3D( ihls, ihls, ihls, ihls, 1, nlay_i )
IF( pszv_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
zsfx_res(ji,jj) = zsfx_res(ji,jj) + pszv_i(ji,jj,jk,jl) * rhoi * z1_dt
pszv_i(ji,jj,jk,jl) = 0._wp
ENDIF
END_3D
ENDDO
ELSE
DO jl = 1, jpl
DO_2D( ihls, ihls, ihls, ihls )
IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
zsfx_res(ji,jj) = zsfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
psv_i (ji,jj,jl) = 0._wp
ENDIF
END_2D
END DO
ENDIF
!
!----------------------------------------
! zap ice energy and send it to the ocean
!----------------------------------------
DO jl = 1, jpl
DO_3D( ihls, ihls, ihls, ihls, 1, nlay_i )
IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
zhfx_res(ji,jj) = zhfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
DO_3D( ihls, ihls, ihls, ihls, 1, nlay_s )
IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_s(ji,jj,jl) <= 0._wp ) THEN
zhfx_res(ji,jj) = zhfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
!--------------------------------------------
! zap ice and snow volume, add water to ocean
!--------------------------------------------
DO jl = 1, jpl
DO_2D( ihls, ihls, ihls, ihls )
IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
zwfx_res(ji,jj) = zwfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
pv_i (ji,jj,jl) = 0._wp
ENDIF
IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
zwfx_res(ji,jj) = zwfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
pv_s (ji,jj,jl) = 0._wp
ENDIF
IF( pv_ip(ji,jj,jl) < 0._wp .OR. pv_il(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN
zwfx_res(ji,jj) = zwfx_res(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt
pv_il (ji,jj,jl) = 0._wp
ENDIF
IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN
zwfx_res(ji,jj) = zwfx_res(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt
pv_ip (ji,jj,jl) = 0._wp
! record residual fluxes
DO_2D( 0, 0, 0, 0 )
wfx_res(ji,jj) = wfx_res(ji,jj) + zwfx_res(ji,jj)
hfx_res(ji,jj) = hfx_res(ji,jj) + zhfx_res(ji,jj)
sfx_res(ji,jj) = sfx_res(ji,jj) + zsfx_res(ji,jj)
END_2D
!
WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
!
END SUBROUTINE ice_var_zapneg
SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i, pszv_i )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_roundoff ***
!!
!! ** Purpose : Remove negative sea ice values arising from roundoff errors
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_il ! melt pond lid volume
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pszv_i ! ice salt content
!!-------------------------------------------------------------------
WHERE( pa_i (1:npti,:) < 0._wp ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
WHERE( pv_i (1:npti,:) < 0._wp ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
WHERE( pv_s (1:npti,:) < 0._wp ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
WHERE( poa_i(1:npti,:) < 0._wp ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
WHERE( pe_i (1:npti,:,:) < 0._wp ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
WHERE( pe_s (1:npti,:,:) < 0._wp ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
WHERE( pa_ip(1:npti,:) < 0._wp ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
WHERE( pv_ip(1:npti,:) < 0._wp ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
WHERE( pv_il(1:npti,:) < 0._wp .AND. pv_il(1:npti,:) > -epsi10 ) pv_il(1:npti,:) = 0._wp ! v_il must be >= 0
IF( nn_icesal == 4 ) THEN
WHERE( pszv_i(1:npti,:,:) < 0._wp ) pszv_i(1:npti,:,:) = 0._wp ! szv_i must be >= 0
ELSE
WHERE( psv_i (1:npti,:) < 0._wp ) psv_i (1:npti,:) = 0._wp ! sv_i must be >= 0
ENDIF
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_brine ***
!! ** Purpose : computes brine volume fraction (%)
!! and salinity of the brine in sea ice
!!
!! ** Method : e = - 0.054 * S (ppt) / T (C)
!!
!! References : Vancoppenolle et al., JGR, 2007
!!-------------------------------------------------------------------
INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: zt1, zt2, zt3, zs_br
!!-------------------------------------------------------------------
!
! brine salinity
zt1 = t_i(ji,jj,jk,jl) - rt0
zt2 = zt1 * zt1
zt3 = zt2 * zt1
IF ( nn_liquidus == 1 ) THEN ; zs_br = - zt1 / rTmlt ! --- Linear liquidus
ELSEIF( nn_liquidus == 2 ) THEN ; zs_br = -18.7_wp * zt1 - 0.519_wp * zt2 - 0.00535_wp * zt3 ! --- 3rd order liquidus, VC19
ELSEIF( nn_liquidus == 3 ) THEN ; zs_br = -17.6_wp * zt1 - 0.389_wp * zt2 - 0.00362_wp * zt3 ! --- Weast 71 liquidus in RJW14
! brine volume fraction
IF( zt1 < - epsi10 ) v_ibr(ji,jj,jl) = v_ibr(ji,jj,jl) + r1_nlay_i * sz_i(ji,jj,jk,jl) / zs_br
DO_2D( 0, 0, 0, 0 )
IF( vt_i(ji,jj) > epsi20 ) THEN
vm_ibr(ji,jj) = SUM( v_ibr(ji,jj,:) * v_i(ji,jj,:) ) / vt_i(ji,jj)
SUBROUTINE ice_var_enthalpy
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_enthalpy ***
!!
!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
!!
!! ** Method : Formula (Bitz and Lipscomb, 1999)
!!-------------------------------------------------------------------
INTEGER :: ji, jk ! dummy loop indices
REAL(wp) :: ztmelts ! local scalar
!!-------------------------------------------------------------------
!
DO jk = 1, nlay_i ! Sea ice energy of melting
DO ji = 1, npti
t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
! (sometimes zdf scheme produces abnormally high temperatures)
e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
& - rcp * ztmelts )
END DO
END DO
DO jk = 1, nlay_s ! Snow energy of melting
DO ji = 1, npti
e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
END DO
END DO
!
END SUBROUTINE ice_var_enthalpy
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
SUBROUTINE ice_var_vremap( ph_old, pts_old, pts_i )
!!-------------------------------------------------------------------
!! *** ROUTINE ice_var_vremap ***
!!
!! ** Purpose :
!! This routine computes new vertical grids in the ice,
!! and consistently redistributes temperatures and salinities
!! Redistribution is made so as to ensure energy/salt conservation
!!
!!
!! ** Method : linear conservative remapping
!!
!! ** Steps : 1) cumulative integrals of old enthalpies/salinities/thicknesses
!! 2) linear remapping on the new layers
!!
!! ------------ cum0(0) ------------- cum1(0)
!! NEW -------------
!! ------------ cum0(1) ==> -------------
!! ... -------------
!! ------------ -------------
!! ------------ cum0(nlay_i+2) ------------- cum1(nlay_i)
!!
!!
!! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005
!!-------------------------------------------------------------------
REAL(wp), DIMENSION(0:nlay_i+1), INTENT(in) :: ph_old, pts_old ! old tickness (m), enthlapy (J.m-2) or salt (m.g/kg)
REAL(wp), DIMENSION(1:nlay_i) , INTENT(inout) :: pts_i ! new enthlapies (J.m-3, remapped) or salt (g/kg)
!
INTEGER :: ji ! dummy loop indices
INTEGER :: jk0, jk1 ! old/new layer indices
!
REAL(wp), DIMENSION(0:nlay_i+2) :: zts_cum0, zh_cum0 ! old cumulative enthlapies/salinities and layers interfaces
REAL(wp), DIMENSION(0:nlay_i) :: zts_cum1, zh_cum1 ! new cumulative enthlapies/salinities and layers interfaces
REAL(wp) :: zhnew ! new layers thicknesses
!!-------------------------------------------------------------------
!-------------------------------------------------------------------------------
! 1) Cumulative integral of old enthalpy/salt * thickness and layers interfaces
!-------------------------------------------------------------------------------
zts_cum0(0) = 0._wp
zh_cum0 (0) = 0._wp
DO jk0 = 1, nlay_i+2
zts_cum0(jk0) = zts_cum0(jk0-1) + pts_old(jk0-1)
zh_cum0 (jk0) = zh_cum0 (jk0-1) + ph_old(jk0-1)
END DO