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
MODULE icedyn_adv_pra
!!======================================================================
!! *** MODULE icedyn_adv_pra ***
!! sea-ice : advection => Prather scheme
!!======================================================================
!! History : ! 2008-03 (M. Vancoppenolle) original code
!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
!!--------------------------------------------------------------------
#if defined key_si3
!!----------------------------------------------------------------------
!! 'key_si3' SI3 sea-ice model
!!----------------------------------------------------------------------
!! ice_dyn_adv_pra : advection of sea ice using Prather scheme
!! adv_x, adv_y : Prather scheme applied in i- and j-direction, resp.
!! adv_pra_init : initialisation of the Prather scheme
!! adv_pra_rst : read/write Prather field in ice restart file, or initialized to zero
!!----------------------------------------------------------------------
USE phycst ! physical constant
USE dom_oce ! ocean domain
USE ice ! sea-ice variables
USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call
USE icevar ! sea-ice: operations
!
USE in_out_manager ! I/O manager
USE iom ! I/O manager library
USE lib_mpp ! MPP library
USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
USE lbclnk ! lateral boundary conditions (or mpp links)
IMPLICIT NONE
PRIVATE
PUBLIC ice_dyn_adv_pra ! called by icedyn_adv
PUBLIC adv_pra_init ! called by icedyn_adv
! Moments for advection
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxage, syage, sxxage, syyage, sxyage ! ice age
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvl , syvl , sxxvl , syyvl , sxyvl ! melt pond lid volume
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/ICE 4.0 , NEMO Consortium (2018)
!! $Id: icedyn_adv_pra.F90 15049 2021-06-23 16:17:30Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
!!----------------------------------------------------------------------
!! ** routine ice_dyn_adv_pra **
!!
!! ** purpose : Computes and adds the advection trend to sea-ice
!!
!! ** method : Uses Prather second order scheme that advects tracers
!! but also their quadratic forms. The method preserves
!! tracer structures by conserving second order moments.
!!
!! Reference: Prather, 1986, JGR, 91, D6. 6671-6681.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
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 thickness
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
!
INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
INTEGER :: icycle ! number of sub-timestep for the advection
REAL(wp) :: zdt, z1_dt ! - -
REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
REAL(wp), DIMENSION(A2D(0)) :: zati1, zati2
REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: ze_s, zes_max
REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp, z0vl
REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es
REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei
!! diagnostics
REAL(wp), DIMENSION(A2D(0)) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
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
!!----------------------------------------------------------------------
!
IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
!
! --- Record max of the surrounding 9-pts (for call Hbig) --- !
! thickness and salinity
WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:)
ELSEWHERE ; zs_i(:,:,:) = 0._wp
END WHERE
CALL icemax3D( ph_i , zhi_max )
CALL icemax3D( ph_s , zhs_max )
CALL icemax3D( ph_ip, zhip_max)
CALL icemax3D( zs_i , zsi_max )
CALL lbc_lnk( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )
!
! enthalpies
DO jk = 1, nlay_i
WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:)
ELSEWHERE ; ze_i(:,:,jk,:) = 0._wp
END WHERE
END DO
DO jk = 1, nlay_s
WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:)
ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp
END WHERE
END DO
CALL icemax4D( ze_i , zei_max )
CALL icemax4D( ze_s , zes_max )
CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp, zes_max, 'T', 1._wp )
!
!
! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
! this should not affect too much the stability
zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rDt_ice * r1_e1u(:,:) )
zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rDt_ice * r1_e2v(:,:) ) )
! non-blocking global communication send zcflnow and receive zcflprv
CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
IF( zcflprv(1) > .5 ) THEN ; icycle = 2
ELSE ; icycle = 1
ENDIF
zdt = rDt_ice / REAL(icycle)
z1_dt = 1._wp / zdt
! --- transport --- !
zudy(:,:) = pu_ice(:,:) * e2u(:,:)
zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
DO jt = 1, icycle
! diagnostics
zdiag_adv_mass(:,:) = SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow
zdiag_adv_salt(:,:) = SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi
zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 )
zati1(:,:) = SUM( pa_i(A2D(0),:), dim=3 )
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
! --- transported fields --- !
DO jl = 1, jpl
zarea(:,:,jl) = e1e2t(:,:)
z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume
z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume
z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area
z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content
z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content
DO jk = 1, nlay_s
z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
END DO
DO jk = 1, nlay_i
z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content
END DO
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction
z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume
IF ( ln_pnd_lids ) THEN
z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:) ! Melt pond lid volume
ENDIF
ENDIF
END DO
!
! !--------------------------------------------!
IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
! !--------------------------------------------!
CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
!
DO jk = 1, nlay_s !--- snow heat content
CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
END DO
DO jk = 1, nlay_i !--- ice heat content
CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
END DO
!
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
IF ( ln_pnd_lids ) THEN
CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
ENDIF
ENDIF
! !--------------------------------------------!
ELSE !== even ice time step: adv_y then adv_x ==!
! !--------------------------------------------!
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
DO jk = 1, nlay_s !--- snow heat content
CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
END DO
DO jk = 1, nlay_i !--- ice heat content
CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
END DO
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
IF ( ln_pnd_lids ) THEN
CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
ENDIF
ENDIF
!
ENDIF
! --- Lateral boundary conditions --- !
! caution: for gradients (sx and sy) the sign changes
CALL lbc_lnk( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp & ! ice volume
& , sxxice, 'T', 1._wp, syyice, 'T', 1._wp, sxyice, 'T', 1._wp &
& , z0snw , 'T', 1._wp, sxsn , 'T', -1._wp, sysn , 'T', -1._wp & ! snw volume
& , sxxsn , 'T', 1._wp, syysn , 'T', 1._wp, sxysn , 'T', 1._wp &
& , z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp & ! ice salinity
& , sxxsal, 'T', 1._wp, syysal, 'T', 1._wp, sxysal, 'T', 1._wp &
& , z0ai , 'T', 1._wp, sxa , 'T', -1._wp, sya , 'T', -1._wp & ! ice concentration
& , sxxa , 'T', 1._wp, syya , 'T', 1._wp, sxya , 'T', 1._wp &
& , z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age
& , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp )
CALL lbc_lnk( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy
& , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp &
& , z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy
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
& , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp )
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
IF( ln_pnd_lids ) THEN
CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction
& , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp &
& , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume
& , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp &
& , z0vl , 'T', 1._wp, sxvl , 'T', -1._wp, syvl , 'T', -1._wp & ! melt pond lid volume
& , sxxvl, 'T', 1._wp, syyvl, 'T', 1._wp, sxyvl, 'T', 1._wp )
ELSE
CALL lbc_lnk( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction
& , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp &
& , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume
& , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp )
ENDIF
ENDIF
! --- Recover the properties from their contents --- !
DO jl = 1, jpl
pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
DO jk = 1, nlay_s
pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
END DO
DO jk = 1, nlay_i
pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
END DO
IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
IF ( ln_pnd_lids ) THEN
pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
ENDIF
ENDIF
END DO
!
! derive open water from ice concentration
zati2(:,:) = SUM( pa_i(A2D(0),:), dim=3 )
pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
& - ( ( zudy(ji,jj) - zudy(ji-1,jj) ) & ! add () for NP repro
& + ( zvdx(ji,jj) - zvdx(ji,jj-1) ) ) * r1_e1e2t(ji,jj) * zdt
END_2D
CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1.0_wp )
!
! --- diagnostics --- !
diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (A2D(0),:) , dim=3 ) * rhoi + SUM( pv_s (A2D(0),:) , dim=3 ) * rhos &
& + SUM( pv_ip(A2D(0),:) , dim=3 ) * rhow + SUM( pv_il(A2D(0),:) , dim=3 ) * rhow &
diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(A2D(0),:) , dim=3 ) * rhoi &
diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(A2D(0),1:nlay_i,:) , dim=4 ), dim=3 ) &
& - SUM(SUM( pe_s(A2D(0),1:nlay_s,:) , dim=4 ), dim=3 ) &
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
& - zdiag_adv_heat(:,:) ) * z1_dt
!
! --- Ensure non-negative fields --- !
! Remove negative values (conservation is ensured)
! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
!
! --- Make sure ice thickness is not too big --- !
! (because ice thickness can be too large where ice concentration is very small)
CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, &
& pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
!
! --- Ensure snow load is not too big --- !
CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
!
END DO
!
IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file
!
END SUBROUTINE ice_dyn_adv_pra
SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , &
& psx, psxx, psy , psyy, psxy )
!!----------------------------------------------------------------------
!! ** routine adv_x **
!!
!! ** purpose : Computes and adds the advection trend to sea-ice
!! variable on x axis
!!----------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! the time step
REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s]
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
!!
INTEGER :: ji, jj, jl, jcat ! dummy loop indices
INTEGER :: jj0 ! dummy loop indices
REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars
REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
REAL(wp) :: zpsm, zps0
REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy
REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - -
REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
!-----------------------------------------------------------------------
! in order to avoid lbc_lnk (communications):
! jj loop must be 1:jpj if adv_x is called first
! and 2:jpj-1 if adv_x is called second
jj0 = NINT(pcrh)
!
jcat = SIZE( ps0 , 3 ) ! size of input arrays
!
DO jl = 1, jcat ! loop on categories
!
! Limitation of moments.
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
DO_2D( 1, 1, jj0, jj0 )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
!
zslpmax = MAX( 0._wp, zps0 )
zs1max = 1.5 * zslpmax
zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
zps0 = zslpmax
zpsx = zs1new * rswitch
zpsxx = zs2new * rswitch
zpsy = zpsy * rswitch
zpsyy = zpsyy * rswitch
zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
! Calculate fluxes and moments between boxes i<-->i+1
! ! Flux from i to i+1 WHEN u GT 0
zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalf1q = zalf1 * zalf1
!
zfm (ji,jj) = zalf * zpsm
zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
zfxx(ji,jj) = zalf * zpsxx * zalfq
zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
zfxy(ji,jj) = zalfq * zpsxy
zfyy(ji,jj) = zalf * zpsyy
! ! Readjust moments remaining in the box.
zpsm = zpsm - zfm(ji,jj)
zps0 = zps0 - zf0(ji,jj)
zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
zpsxx = zalf1 * zalf1q * zpsxx
zpsy = zpsy - zfy (ji,jj)
zpsyy = zpsyy - zfyy(ji,jj)
zpsxy = zalf1q * zpsxy
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
!
END_2D
DO_2D( 1, 0, jj0, jj0 )
! ! Flux from i+1 to i when u LT 0.
zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
zalg (ji,jj) = zalf
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalg1 (ji,jj) = zalf1
zalf1q = zalf1 * zalf1
zalg1q(ji,jj) = zalf1q
!
zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
& - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
zfxx (ji,jj) = zfxx(ji,jj) + zalf * ( psxx(ji+1,jj,jl) * zalfq )
zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
END_2D
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
! ! Readjust moments remaining in the box.
zbt = zbet(ji-1,jj)
zbt1 = 1.0 - zbet(ji-1,jj)
!
zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
zpsxy = zalg1q(ji-1,jj) * zpsxy
! Put the temporary moments into appropriate neighboring boxes.
! ! Flux from i to i+1 IF u GT 0.
zbt = zbet(ji-1,jj)
zbt1 = 1.0 - zbet(ji-1,jj)
zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji-1,jj) )
zalf = zbt * zfm(ji-1,jj) / zpsm
zalf1 = 1.0 - zalf
ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
!
zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji-1,jj) )
zpsx = zbt1 * zpsx + zbt * ( ( zalf * zfx(ji-1,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
zpsxx = zbt1 * zpsxx + zbt * ( ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx ) &
& + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) )
zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy ) &
& + 3.0 * ( -zalf1 * zfy(ji-1,jj) + zalf * zpsy ) )
zpsy = zbt1 * zpsy + zbt * ( zpsy + zfy (ji-1,jj) )
zpsyy = zbt1 * zpsyy + zbt * ( zpsyy + zfyy(ji-1,jj) )
! ! Flux from i+1 to i IF u LT 0.
zbt = zbet(ji,jj)
zbt1 = 1.0 - zbet(ji,jj)
zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji ,jj) )
zalf = zbt1 * zfm(ji ,jj) / zpsm
zalf1 = 1.0 - zalf
ztemp = -zalf * zps0 + zalf1 * zf0(ji ,jj)
!
zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji ,jj) )
zpsx = zbt * zpsx + zbt1 * ( ( zalf * zfx(ji ,jj) + zalf1 * zpsx ) + 3.0 * ztemp )
zpsxx = zbt * zpsxx + zbt1 * ( ( zalf * zalf * zfxx(ji ,jj) + zalf1 * zalf1 * zpsxx ) &
& + 5.0 * ( zalf * zalf1 * ( -zpsx + zfx(ji ,jj) ) + ( zalf1 - zalf ) * ztemp ) )
zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji ,jj) + zalf1 * zpsxy ) &
& + 3.0 * ( zalf1 * zfy(ji ,jj) - zalf * zpsy ) )
zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji ,jj) )
zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji ,jj) )
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
!
END_2D
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
!
END DO
!
END SUBROUTINE adv_x
SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , &
& psx, psxx, psy , psyy, psxy )
!!---------------------------------------------------------------------
!! ** routine adv_y **
!!
!! ** purpose : Computes and adds the advection trend to sea-ice
!! variable on y axis
!!---------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! time step
REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s]
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
!!
INTEGER :: ji, jj, jl, jcat ! dummy loop indices
INTEGER :: ji0 ! dummy loop indices
REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars
REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
REAL(wp) :: zpsm, zps0
REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy
REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - -
REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
!---------------------------------------------------------------------
! in order to avoid lbc_lnk (communications):
! ji loop must be 1:jpi if adv_y is called first
! and 2:jpi-1 if adv_y is called second
ji0 = NINT(pcrh)
!
jcat = SIZE( ps0 , 3 ) ! size of input arrays
!
DO jl = 1, jcat ! loop on categories
!
! Limitation of moments.
DO_2D( ji0, ji0, 1, 1 )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
!
! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise)
zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
!
zslpmax = MAX( 0._wp, zps0 )
zs1max = 1.5 * zslpmax
zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) )
zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new ) - zslpmax, zpsyy ) )
rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
!
zps0 = zslpmax
zpsx = zpsx * rswitch
zpsxx = zpsxx * rswitch
zpsy = zs1new * rswitch
zpsyy = zs2new * rswitch
zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
! Calculate fluxes and moments between boxes j<-->j+1
! ! Flux from j to j+1 WHEN v GT 0
zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalf1q = zalf1 * zalf1
!
zfm (ji,jj) = zalf * zpsm
zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1 - zalf) * zpsyy ) )
zfy (ji,jj) = zalfq * ( zpsy + 3.0 * zalf1 * zpsyy )
zfyy(ji,jj) = zalf * zpsyy * zalfq
zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
zfxy(ji,jj) = zalfq * zpsxy
zfxx(ji,jj) = zalf * zpsxx
!
! ! Readjust moments remaining in the box.
zpsm = zpsm - zfm(ji,jj)
zps0 = zps0 - zf0(ji,jj)
zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy )
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
zpsx = zpsx - zfx(ji,jj)
zpsxx = zpsxx - zfxx(ji,jj)
zpsxy = zalf1q * zpsxy
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
END_2D
!
DO_2D( ji0, ji0, 1, 0 )
! ! Flux from j+1 to j when v LT 0.
zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)
zalg (ji,jj) = zalf
zalfq = zalf * zalf
zalf1 = 1.0 - zalf
zalg1 (ji,jj) = zalf1
zalf1q = zalf1 * zalf1
zalg1q(ji,jj) = zalf1q
!
zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl)
zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) &
& - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
zfyy (ji,jj) = zfyy(ji,jj) + zalf * ( psyy(ji,jj+1,jl) * zalfq )
zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl)
zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl)
END_2D
DO_2D( ji0, ji0, 0, 0 )
!
zpsm = psm (ji,jj,jl) ! optimization
zps0 = ps0 (ji,jj,jl)
zpsx = psx (ji,jj,jl)
zpsxx = psxx(ji,jj,jl)
zpsy = psy (ji,jj,jl)
zpsyy = psyy(ji,jj,jl)
zpsxy = psxy(ji,jj,jl)
! ! Readjust moments remaining in the box.
zbt = zbet(ji,jj-1)
zbt1 = ( 1.0 - zbet(ji,jj-1) )
zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
zpsxy = zalg1q(ji,jj-1) * zpsxy
! Put the temporary moments into appropriate neighboring boxes.
! ! Flux from j to j+1 IF v GT 0.
zbt = zbet(ji,jj-1)
zbt1 = 1.0 - zbet(ji,jj-1)
zpsm = zbt1 * zpsm + zbt * ( zpsm + zfm(ji,jj-1) )
zalf = zbt * zfm(ji,jj-1) / zpsm
ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
zps0 = zbt1 * zps0 + zbt * ( zps0 + zf0(ji,jj-1) )
zpsy = zbt1 * zpsy + zbt * ( ( zalf * zfy(ji,jj-1) + zalf1 * zpsy ) + 3.0 * ztemp )
zpsyy = zbt1 * zpsyy + zbt * ( ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy ) &
& + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )
zpsxy = zbt1 * zpsxy + zbt * ( ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy ) &
& + 3.0 * ( -zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )
zpsx = zbt1 * zpsx + zbt * ( zpsx + zfx (ji,jj-1) )
zpsxx = zbt1 * zpsxx + zbt * ( zpsxx + zfxx(ji,jj-1) )
! ! Flux from j+1 to j IF v LT 0.
zbt = zbet(ji,jj)
zbt1 = 1.0 - zbet(ji,jj)
zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
zalf = zbt1 * zfm(ji,jj) / zpsm
ztemp = -zalf * zps0 + zalf1 * zf0(ji,jj)
zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj ) )
zpsy = zbt * zpsy + zbt1 * ( ( zalf * zfy(ji,jj ) + zalf1 * zpsy ) + 3.0 * ztemp )
zpsyy = zbt * zpsyy + zbt1 * ( ( zalf * zalf * zfyy(ji,jj ) + zalf1 * zalf1 * zpsyy ) &
& + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj ) ) + ( zalf1 - zalf ) * ztemp ) )
zpsxy = zbt * zpsxy + zbt1 * ( ( zalf * zfxy(ji,jj ) + zalf1 * zpsxy ) &
& + 3.0 * ( zalf1 * zfx(ji,jj ) - zalf * zpsx ) )
zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj ) )
zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj ) )
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
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
!
psm (ji,jj,jl) = zpsm ! optimization
ps0 (ji,jj,jl) = zps0
psx (ji,jj,jl) = zpsx
psxx(ji,jj,jl) = zpsxx
psy (ji,jj,jl) = zpsy
psyy(ji,jj,jl) = zpsyy
psxy(ji,jj,jl) = zpsxy
END_2D
!
END DO
!
END SUBROUTINE adv_y
SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, &
& pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
!!-------------------------------------------------------------------
!! *** ROUTINE Hbig ***
!!
!! ** Purpose : Thickness correction in case advection scheme creates
!! abnormally tick ice or snow
!!
!! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
!! (before advection) and reduce it by adapting ice concentration
!! 2- check whether snow thickness is larger than the surrounding 9-points
!! (before advection) and reduce it by sending the excess in the ocean
!!
!! ** input : Max thickness of the surrounding 9-points
!!-------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max
REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i
!
INTEGER :: ji, jj, jk, jl ! dummy loop indices
REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra
!!-------------------------------------------------------------------
!
z1_dt = 1._wp / pdt
!
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
!
! ! -- check h_ip -- !
! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
ENDIF
ENDIF
!
! ! -- check h_i -- !
! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)
ENDIF
!
! ! -- check h_s -- !
! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
!
wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt
hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
!
pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
ENDIF
!
! ! -- check s_i -- !
! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = psi_max(ji,jj,jl) / zsi
sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt
psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra
ENDIF
!
ENDIF
END_2D
END DO
!
! ! -- check e_i/v_i -- !
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i )
IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = pei_max(ji,jj,jk,jl) / zei
hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra
ENDIF
ENDIF
END_3D
END DO
! ! -- check e_s/v_s -- !
DO jl = 1, jpl
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s )
IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
zfra = pes_max(ji,jj,jk,jl) / zes
hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra
ENDIF
ENDIF
END_3D
END DO
!
END SUBROUTINE Hbig
SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
!!-------------------------------------------------------------------
!! *** ROUTINE Hsnow ***
!!
!! ** Purpose : 1- Check snow load after advection
!! 2- Correct pond concentration to avoid a_ip > a_i
!!
!! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface
!! then put the snow excess in the ocean
!!
!! ** Notes : This correction is crucial because of the call to routine icecor afterwards
!! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
!! make the snow very thick (if concentration decreases drastically)
!! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
!!-------------------------------------------------------------------
REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip
REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
!
INTEGER :: ji, jj, jl ! dummy loop indices
REAL(wp) :: z1_dt, zvs_excess, zfra
!!-------------------------------------------------------------------
!
z1_dt = 1._wp / pdt
!
! -- check snow load -- !
DO jl = 1, jpl
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
!
zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rho0-rhoi) * r1_rhos )
!
IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface
! put snow excess in the ocean
zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
! correct snow volume and heat content
pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
ENDIF
!
ENDIF
END_2D
END DO
!
!-- correct pond concentration to avoid a_ip > a_i -- !
WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
!
END SUBROUTINE Hsnow
SUBROUTINE adv_pra_init
!!-------------------------------------------------------------------
!! *** ROUTINE adv_pra_init ***
!!
!! ** Purpose : allocate and initialize arrays for Prather advection
!!-------------------------------------------------------------------
INTEGER :: ierr
!!-------------------------------------------------------------------
!
! !* allocate prather fields
ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , &
& sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , &
& sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , &
& sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) , &
& sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) , &
& sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , &
& sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , &
& sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) , &
!
& sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
& syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl) , &
!
& sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
& syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl) , &
& STAT = ierr )
!
CALL mpp_sum( 'icedyn_adv_pra', ierr )
IF( ierr /= 0 ) CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
!
CALL adv_pra_rst( 'READ' ) !* read or initialize all required files
!
END SUBROUTINE adv_pra_init
SUBROUTINE adv_pra_rst( cdrw, kt )
!!---------------------------------------------------------------------
!! *** ROUTINE adv_pra_rst ***
!!
!! ** Purpose : Read or write file in restart file
!!
!! ** Method : use of IOM library
!!----------------------------------------------------------------------
CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag
INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
!
INTEGER :: jk, jl ! dummy loop indices
INTEGER :: iter ! local integer
INTEGER :: id1 ! local integer
CHARACTER(len=25) :: znam
CHARACTER(len=2) :: zchar, zchar1
REAL(wp), DIMENSION(jpi,jpj,jpl) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
! !==========================!
IF( TRIM(cdrw) == 'READ' ) THEN !== Read or initialize ==!
! !==========================!
!
IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0
ELSE ; id1 = 0 ! no restart: id1=0
ENDIF
!
IF( id1 > 0 ) THEN !** Read the restart file **!
!
! ! ice thickness
CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice )
CALL iom_get( numrir, jpdom_auto, 'syyice', syyice )
CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice )
! ! snow thickness
CALL iom_get( numrir, jpdom_auto, 'sxsn' , sxsn , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sysn' , sysn , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn )
CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn )
CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn )
! ! ice concentration
CALL iom_get( numrir, jpdom_auto, 'sxa' , sxa , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sya' , sya , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxa' , sxxa )
CALL iom_get( numrir, jpdom_auto, 'syya' , syya )
CALL iom_get( numrir, jpdom_auto, 'sxya' , sxya )
! ! ice salinity
CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal )
CALL iom_get( numrir, jpdom_auto, 'syysal', syysal )
CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal )
! ! ice age
CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp )
CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage )
CALL iom_get( numrir, jpdom_auto, 'syyage', syyage )
CALL iom_get( numrir, jpdom_auto, 'sxyage', sxyage )
! ! snow layers heat content
DO jk = 1, nlay_s
WRITE(zchar1,'(I2.2)') jk
znam = 'sxc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)
znam = 'syc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp ) ; syc0 (:,:,jk,:) = z3d(:,:,:)
znam = 'sxxc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:)
znam = 'syyc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:)
znam = 'sxyc0'//'_l'//zchar1
CALL iom_get( numrir, jpdom_auto, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:)
END DO
! ! ice layers heat content
DO jk = 1, nlay_i