Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icedyn_adv_pra.F90

Last change on this file was 14026, checked in by clem, 4 years ago
4.0-HEAD: fix bugs and defects related to tickets #2573 #2575 #2576 #2578. Sette passed and those fixes are now in the trunk, so unless there is a tricky trick somewhere, everything should be fine.
Property svn:keywords set to `Id`
File size: 71.2 KB

Line
1	MODULE icedyn_adv_pra
2	!!======================================================================
3	!! * MODULE icedyn_adv_pra *
4	!! sea-ice : advection => Prather scheme
5	!!======================================================================
6	!! History : ! 2008-03 (M. Vancoppenolle) original code
7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
8	!!--------------------------------------------------------------------
9	#if defined key_si3
10	!!----------------------------------------------------------------------
11	!! 'key_si3' SI3 sea-ice model
12	!!----------------------------------------------------------------------
13	!! ice_dyn_adv_pra : advection of sea ice using Prather scheme
14	!! adv_x, adv_y : Prather scheme applied in i- and j-direction, resp.
15	!! adv_pra_init : initialisation of the Prather scheme
16	!! adv_pra_rst : read/write Prather field in ice restart file, or initialized to zero
17	!!----------------------------------------------------------------------
18	USE phycst ! physical constant
19	USE dom_oce ! ocean domain
20	USE ice ! sea-ice variables
21	USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call
22	USE icevar ! sea-ice: operations
23	!
24	USE in_out_manager ! I/O manager
25	USE iom ! I/O manager library
26	USE lib_mpp ! MPP library
27	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
28	USE lbclnk ! lateral boundary conditions (or mpp links)
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC ice_dyn_adv_pra ! called by icedyn_adv
34	PUBLIC adv_pra_init ! called by icedyn_adv
35
36	! Moments for advection
37	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness
38	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness
39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration
40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity
41	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxage, syage, sxxage, syyage, sxyage ! ice age
42	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content
43	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content
44	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction
45	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume
46	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvl , syvl , sxxvl , syyvl , sxyvl ! melt pond lid volume
47
48	!! * Substitutions
49	# include "vectopt_loop_substitute.h90"
50	!!----------------------------------------------------------------------
51	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
52	!! $Id$
53	!! Software governed by the CeCILL license (see ./LICENSE)
54	!!----------------------------------------------------------------------
55	CONTAINS
56
57	SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, &
58	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i )
59	!!----------------------------------------------------------------------
60	!! routine ice_dyn_adv_pra
61	!!
62	!! ** purpose : Computes and adds the advection trend to sea-ice
63	!!
64	!! ** method : Uses Prather second order scheme that advects tracers
65	!! but also their quadratic forms. The method preserves
66	!! tracer structures by conserving second order moments.
67	!!
68	!! Reference: Prather, 1986, JGR, 91, D6. 6671-6681.
69	!!----------------------------------------------------------------------
70	INTEGER , INTENT(in ) :: kt ! time step
71	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
72	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
73	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness
74	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness
75	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness
76	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
77	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
78	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
79	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
80	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
81	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
82	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
83	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
84	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_il ! melt pond lid thickness
85	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
86	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
87	!
88	INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
89	INTEGER :: icycle ! number of sub-timestep for the advection
90	REAL(wp) :: zdt, z1_dt ! - -
91	REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
92	REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
93	REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx
94	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max, zs_i, zsi_max
95	REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: ze_i, zei_max
96	REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: ze_s, zes_max
97	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea
98	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi
99	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp, z0vl
100	REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es
101	REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei
102	!! diagnostics
103	REAL(wp), DIMENSION(jpi,jpj) :: zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
104	!!----------------------------------------------------------------------
105	!
106	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
107	!
108	! --- Record max of the surrounding 9-pts (for call Hbig) --- !
109	! thickness and salinity
110	WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:)
111	ELSEWHERE ; zs_i(:,:,:) = 0._wp
112	END WHERE
113	CALL icemax3D( ph_i , zhi_max )
114	CALL icemax3D( ph_s , zhs_max )
115	CALL icemax3D( ph_ip, zhip_max)
116	CALL icemax3D( zs_i , zsi_max )
117	CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. )
118
119	! enthalpies
120	DO jk = 1, nlay_i
121	WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:)
122	ELSEWHERE ; ze_i(:,:,jk,:) = 0._wp
123	END WHERE
124	END DO
125	DO jk = 1, nlay_s
126	WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:)
127	ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp
128	END WHERE
129	END DO
130	CALL icemax4D( ze_i , zei_max )
131	CALL icemax4D( ze_s , zes_max )
132	CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. )
133	CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. )
134	!
135	!
136	! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
137	! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
138	! this should not affect too much the stability
139	zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
140	zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
141
142	! non-blocking global communication send zcflnow and receive zcflprv
143	CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
144
145	IF( zcflprv(1) > .5 ) THEN ; icycle = 2
146	ELSE ; icycle = 1
147	ENDIF
148	zdt = rdt_ice / REAL(icycle)
149	z1_dt = 1._wp / zdt
150
151	! --- transport --- !
152	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
153	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
154
155	DO jt = 1, icycle
156
157	! diagnostics
158	zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
159	& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow
160	zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi
161	zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
162	& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
163
164	! record at_i before advection (for open water)
165	zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
166
167	! --- transported fields --- !
168	DO jl = 1, jpl
169	zarea(:,:,jl) = e1e2t(:,:)
170	z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume
171	z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume
172	z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area
173	z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content
174	z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content
175	DO jk = 1, nlay_s
176	z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
177	END DO
178	DO jk = 1, nlay_i
179	z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content
180	END DO
181	IF ( ln_pnd_LEV ) THEN
182	z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction
183	z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume
184	IF ( ln_pnd_lids ) THEN
185	z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:) ! Melt pond lid volume
186	ENDIF
187	ENDIF
188	END DO
189	!
190	! !--------------------------------------------!
191	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
192	! !--------------------------------------------!
193	CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
194	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
195	CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
196	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
197	CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
198	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
199	CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
200	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
201	CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
202	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
203	!
204	DO jk = 1, nlay_s !--- snow heat content
205	CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
206	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
207	CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
208	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
209	END DO
210	DO jk = 1, nlay_i !--- ice heat content
211	CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
212	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
213	CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
214	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
215	END DO
216	!
217	IF ( ln_pnd_LEV ) THEN
218	CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
219	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
220	CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
221	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
222	IF ( ln_pnd_lids ) THEN
223	CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
224	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
225	ENDIF
226	ENDIF
227	! !--------------------------------------------!
228	ELSE !== even ice time step: adv_y then adv_x ==!
229	! !--------------------------------------------!
230	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
231	CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
232	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
233	CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
234	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
235	CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
236	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
237	CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
238	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
239	CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
240	DO jk = 1, nlay_s !--- snow heat content
241	CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
242	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
243	CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
244	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
245	END DO
246	DO jk = 1, nlay_i !--- ice heat content
247	CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
248	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
249	CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
250	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
251	END DO
252	IF ( ln_pnd_LEV ) THEN
253	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
254	CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
255	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
256	CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
257	IF ( ln_pnd_lids ) THEN
258	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume
259	CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )
260	ENDIF
261	ENDIF
262	!
263	ENDIF
264
265	! --- Lateral boundary conditions --- !
266	! caution: for gradients (sx and sy) the sign changes
267	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp & ! ice volume
268	& , sxxice, 'T', 1._wp, syyice, 'T', 1._wp, sxyice, 'T', 1._wp &
269	& , z0snw , 'T', 1._wp, sxsn , 'T', -1._wp, sysn , 'T', -1._wp & ! snw volume
270	& , sxxsn , 'T', 1._wp, syysn , 'T', 1._wp, sxysn , 'T', 1._wp )
271	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp & ! ice salinity
272	& , sxxsal, 'T', 1._wp, syysal, 'T', 1._wp, sxysal, 'T', 1._wp &
273	& , z0ai , 'T', 1._wp, sxa , 'T', -1._wp, sya , 'T', -1._wp & ! ice concentration
274	& , sxxa , 'T', 1._wp, syya , 'T', 1._wp, sxya , 'T', 1._wp )
275	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp & ! ice age
276	& , sxxage, 'T', 1._wp, syyage, 'T', 1._wp, sxyage, 'T', 1._wp )
277	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es , 'T', 1._wp, sxc0 , 'T', -1._wp, syc0 , 'T', -1._wp & ! snw enthalpy
278	& , sxxc0 , 'T', 1._wp, syyc0 , 'T', 1._wp, sxyc0 , 'T', 1._wp )
279	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy
280	& , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp )
281	IF ( ln_pnd_LEV ) THEN
282	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction
283	& , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp &
284	& , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp & ! melt pond volume
285	& , sxxvp, 'T', 1._wp, syyvp, 'T', 1._wp, sxyvp, 'T', 1._wp )
286	IF ( ln_pnd_lids ) THEN
287	CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp & ! melt pond lid volume
288	& , sxxvl,'T', 1._wp, syyvl,'T', 1._wp, sxyvl,'T', 1._wp )
289	ENDIF
290	ENDIF
291
292	! --- Recover the properties from their contents --- !
293	DO jl = 1, jpl
294	pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
295	pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
296	psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
297	poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
298	pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
299	DO jk = 1, nlay_s
300	pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
301	END DO
302	DO jk = 1, nlay_i
303	pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
304	END DO
305	IF ( ln_pnd_LEV ) THEN
306	pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
307	pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
308	IF ( ln_pnd_lids ) THEN
309	pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
310	ENDIF
311	ENDIF
312	END DO
313	!
314	! derive open water from ice concentration
315	zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
316	DO jj = 2, jpjm1
317	DO ji = fs_2, fs_jpim1
318	pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
319	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
320	END DO
321	END DO
322	CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. )
323	!
324	! --- diagnostics --- !
325	diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos &
326	& + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow &
327	& - zdiag_adv_mass(:,:) ) * z1_dt
328	diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
329	& - zdiag_adv_salt(:,:) ) * z1_dt
330	diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
331	& - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
332	& - zdiag_adv_heat(:,:) ) * z1_dt
333	!
334	! --- Ensure non-negative fields --- !
335	! Remove negative values (conservation is ensured)
336	! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
337	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 )
338	!
339	! --- Make sure ice thickness is not too big --- !
340	! (because ice thickness can be too large where ice concentration is very small)
341	CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, &
342	& pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
343	!
344	! --- Ensure snow load is not too big --- !
345	CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
346	!
347	END DO
348	!
349	IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file
350	!
351	END SUBROUTINE ice_dyn_adv_pra
352
353
354	SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , &
355	& psx, psxx, psy , psyy, psxy )
356	!!----------------------------------------------------------------------
357	!! routine adv_x
358	!!
359	!! ** purpose : Computes and adds the advection trend to sea-ice
360	!! variable on x axis
361	!!----------------------------------------------------------------------
362	REAL(wp) , INTENT(in ) :: pdt ! the time step
363	REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
364	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s]
365	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
366	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
367	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
368	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
369	!!
370	INTEGER :: ji, jj, jl, jcat ! dummy loop indices
371	INTEGER :: jjmin, jjmax ! dummy loop indices
372	REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars
373	REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
374	REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
375	REAL(wp) :: zpsm, zps0
376	REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy
377	REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace
378	REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - -
379	REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
380	!-----------------------------------------------------------------------
381	! in order to avoid lbc_lnk (communications):
382	! jj loop must be 1:jpj if adv_x is called first
383	! and 2:jpj-1 if adv_x is called second
384	jjmin = 2 - NINT(pcrh) ! 1 or 2
385	jjmax = jpjm1 + NINT(pcrh) ! jpj or jpj-1
386	!
387	jcat = SIZE( ps0 , 3 ) ! size of input arrays
388	!
389	DO jl = 1, jcat ! loop on categories
390	!
391	! Limitation of moments.
392	DO jj = jjmin, jjmax
393
394	DO ji = 1, jpi
395
396	zpsm = psm (ji,jj,jl) ! optimization
397	zps0 = ps0 (ji,jj,jl)
398	zpsx = psx (ji,jj,jl)
399	zpsxx = psxx(ji,jj,jl)
400	zpsy = psy (ji,jj,jl)
401	zpsyy = psyy(ji,jj,jl)
402	zpsxy = psxy(ji,jj,jl)
403
404	! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
405	zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
406	!
407	zslpmax = MAX( 0._wp, zps0 )
408	zs1max = 1.5 * zslpmax
409	zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) )
410	zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
411	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
412
413	zps0 = zslpmax
414	zpsx = zs1new * rswitch
415	zpsxx = zs2new * rswitch
416	zpsy = zpsy * rswitch
417	zpsyy = zpsyy * rswitch
418	zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
419
420	! Calculate fluxes and moments between boxes i<-->i+1
421	! ! Flux from i to i+1 WHEN u GT 0
422	zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
423	zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
424	zalfq = zalf * zalf
425	zalf1 = 1.0 - zalf
426	zalf1q = zalf1 * zalf1
427	!
428	zfm (ji,jj) = zalf * zpsm
429	zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
430	zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx )
431	zfxx(ji,jj) = zalf * zpsxx * zalfq
432	zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy )
433	zfxy(ji,jj) = zalfq * zpsxy
434	zfyy(ji,jj) = zalf * zpsyy
435
436	! ! Readjust moments remaining in the box.
437	zpsm = zpsm - zfm(ji,jj)
438	zps0 = zps0 - zf0(ji,jj)
439	zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
440	zpsxx = zalf1 * zalf1q * zpsxx
441	zpsy = zpsy - zfy (ji,jj)
442	zpsyy = zpsyy - zfyy(ji,jj)
443	zpsxy = zalf1q * zpsxy
444	!
445	psm (ji,jj,jl) = zpsm ! optimization
446	ps0 (ji,jj,jl) = zps0
447	psx (ji,jj,jl) = zpsx
448	psxx(ji,jj,jl) = zpsxx
449	psy (ji,jj,jl) = zpsy
450	psyy(ji,jj,jl) = zpsyy
451	psxy(ji,jj,jl) = zpsxy
452	!
453	END DO
454
455	DO ji = 1, fs_jpim1
456	! ! Flux from i+1 to i when u LT 0.
457	zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
458	zalg (ji,jj) = zalf
459	zalfq = zalf * zalf
460	zalf1 = 1.0 - zalf
461	zalg1 (ji,jj) = zalf1
462	zalf1q = zalf1 * zalf1
463	zalg1q(ji,jj) = zalf1q
464	!
465	zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
466	zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
467	& - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
468	zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
469	zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq
470	zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
471	zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
472	zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
473	END DO
474
475	DO ji = fs_2, fs_jpim1
476	!
477	zpsm = psm (ji,jj,jl) ! optimization
478	zps0 = ps0 (ji,jj,jl)
479	zpsx = psx (ji,jj,jl)
480	zpsxx = psxx(ji,jj,jl)
481	zpsy = psy (ji,jj,jl)
482	zpsyy = psyy(ji,jj,jl)
483	zpsxy = psxy(ji,jj,jl)
484	! ! Readjust moments remaining in the box.
485	zbt = zbet(ji-1,jj)
486	zbt1 = 1.0 - zbet(ji-1,jj)
487	!
488	zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
489	zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
490	zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
491	zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
492	zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) )
493	zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
494	zpsxy = zalg1q(ji-1,jj) * zpsxy
495
496	! Put the temporary moments into appropriate neighboring boxes.
497	! ! Flux from i to i+1 IF u GT 0.
498	zbt = zbet(ji-1,jj)
499	zbt1 = 1.0 - zbet(ji-1,jj)
500	zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
501	zalf = zbt * zfm(ji-1,jj) / zpsm
502	zalf1 = 1.0 - zalf
503	ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
504	!
505	zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
506	zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx
507	zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx &
508	& + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
509	& + zbt1 * zpsxx
510	zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy &
511	& + 3.0 * (- zalf1zfy(ji-1,jj) + zalf zpsy ) ) &
512	& + zbt1 * zpsxy
513	zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy
514	zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
515
516	! ! Flux from i+1 to i IF u LT 0.
517	zbt = zbet(ji,jj)
518	zbt1 = 1.0 - zbet(ji,jj)
519	zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
520	zalf = zbt1 * zfm(ji,jj) / zpsm
521	zalf1 = 1.0 - zalf
522	ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
523	!
524	zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
525	zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp )
526	zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx &
527	& + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) &
528	& + ( zalf1 - zalf ) * ztemp ) )
529	zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy &
530	& + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
531	zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) )
532	zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
533	!
534	psm (ji,jj,jl) = zpsm ! optimization
535	ps0 (ji,jj,jl) = zps0
536	psx (ji,jj,jl) = zpsx
537	psxx(ji,jj,jl) = zpsxx
538	psy (ji,jj,jl) = zpsy
539	psyy(ji,jj,jl) = zpsyy
540	psxy(ji,jj,jl) = zpsxy
541	END DO
542
543	END DO
544
545	END DO
546	!
547	END SUBROUTINE adv_x
548
549
550	SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , &
551	& psx, psxx, psy , psyy, psxy )
552	!!---------------------------------------------------------------------
553	!! routine adv_y
554	!!
555	!! ** purpose : Computes and adds the advection trend to sea-ice
556	!! variable on y axis
557	!!---------------------------------------------------------------------
558	REAL(wp) , INTENT(in ) :: pdt ! time step
559	REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
560	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s]
561	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
562	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
563	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
564	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
565	!!
566	INTEGER :: ji, jj, jl, jcat ! dummy loop indices
567	INTEGER :: jimin, jimax ! dummy loop indices
568	REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars
569	REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
570	REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
571	REAL(wp) :: zpsm, zps0
572	REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy
573	REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace
574	REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - -
575	REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
576	!---------------------------------------------------------------------
577	! in order to avoid lbc_lnk (communications):
578	! ji loop must be 1:jpi if adv_y is called first
579	! and 2:jpi-1 if adv_y is called second
580	jimin = 2 - NINT(pcrh) ! 1 or 2
581	jimax = jpim1 + NINT(pcrh) ! jpi or jpi-1
582	!
583	jcat = SIZE( ps0 , 3 ) ! size of input arrays
584	!
585	DO jl = 1, jcat ! loop on categories
586	!
587	! Limitation of moments.
588	DO jj = 1, jpj
589	DO ji = jimin, jimax
590	!
591	zpsm = psm (ji,jj,jl) ! optimization
592	zps0 = ps0 (ji,jj,jl)
593	zpsx = psx (ji,jj,jl)
594	zpsxx = psxx(ji,jj,jl)
595	zpsy = psy (ji,jj,jl)
596	zpsyy = psyy(ji,jj,jl)
597	zpsxy = psxy(ji,jj,jl)
598	!
599	! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise)
600	zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
601	!
602	zslpmax = MAX( 0._wp, zps0 )
603	zs1max = 1.5 * zslpmax
604	zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) )
605	zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
606	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
607	!
608	zps0 = zslpmax
609	zpsx = zpsx * rswitch
610	zpsxx = zpsxx * rswitch
611	zpsy = zs1new * rswitch
612	zpsyy = zs2new * rswitch
613	zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
614
615	! Calculate fluxes and moments between boxes j<-->j+1
616	! ! Flux from j to j+1 WHEN v GT 0
617	zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
618	zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm
619	zalfq = zalf * zalf
620	zalf1 = 1.0 - zalf
621	zalf1q = zalf1 * zalf1
622	!
623	zfm (ji,jj) = zalf * zpsm
624	zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) )
625	zfy (ji,jj) = zalfq ( zpsy + 3.0zalf1*zpsyy )
626	zfyy(ji,jj) = zalf * zalfq * zpsyy
627	zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy )
628	zfxy(ji,jj) = zalfq * zpsxy
629	zfxx(ji,jj) = zalf * zpsxx
630	!
631	! ! Readjust moments remaining in the box.
632	zpsm = zpsm - zfm(ji,jj)
633	zps0 = zps0 - zf0(ji,jj)
634	zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy )
635	zpsyy = zalf1 * zalf1q * zpsyy
636	zpsx = zpsx - zfx(ji,jj)
637	zpsxx = zpsxx - zfxx(ji,jj)
638	zpsxy = zalf1q * zpsxy
639	!
640	psm (ji,jj,jl) = zpsm ! optimization
641	ps0 (ji,jj,jl) = zps0
642	psx (ji,jj,jl) = zpsx
643	psxx(ji,jj,jl) = zpsxx
644	psy (ji,jj,jl) = zpsy
645	psyy(ji,jj,jl) = zpsyy
646	psxy(ji,jj,jl) = zpsxy
647	END DO
648	END DO
649	!
650	DO jj = 1, jpjm1
651	DO ji = jimin, jimax
652	! ! Flux from j+1 to j when v LT 0.
653	zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)
654	zalg (ji,jj) = zalf
655	zalfq = zalf * zalf
656	zalf1 = 1.0 - zalf
657	zalg1 (ji,jj) = zalf1
658	zalf1q = zalf1 * zalf1
659	zalg1q(ji,jj) = zalf1q
660	!
661	zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl)
662	zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) &
663	& - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
664	zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
665	zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq
666	zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
667	zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl)
668	zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl)
669	END DO
670	END DO
671
672	DO jj = 2, jpjm1
673	DO ji = jimin, jimax
674	! ! Readjust moments remaining in the box.
675	zbt = zbet(ji,jj-1)
676	zbt1 = ( 1.0 - zbet(ji,jj-1) )
677	!
678	zpsm = psm (ji,jj,jl) ! optimization
679	zps0 = ps0 (ji,jj,jl)
680	zpsx = psx (ji,jj,jl)
681	zpsxx = psxx(ji,jj,jl)
682	zpsy = psy (ji,jj,jl)
683	zpsyy = psyy(ji,jj,jl)
684	zpsxy = psxy(ji,jj,jl)
685	!
686	zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
687	zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
688	zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
689	zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
690	zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) )
691	zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
692	zpsxy = zalg1q(ji,jj-1) * zpsxy
693
694	! Put the temporary moments into appropriate neighboring boxes.
695	! ! Flux from j to j+1 IF v GT 0.
696	zbt = zbet(ji,jj-1)
697	zbt1 = 1.0 - zbet(ji,jj-1)
698	zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm
699	zalf = zbt * zfm(ji,jj-1) / zpsm
700	zalf1 = 1.0 - zalf
701	ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
702	!
703	zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
704	zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) &
705	& + zbt1 * zpsy
706	zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy &
707	& + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &
708	& + zbt1 * zpsyy
709	zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy &
710	& + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) &
711	& + zbt1 * zpsxy
712	zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx
713	zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
714
715	! ! Flux from j+1 to j IF v LT 0.
716	zbt = zbet(ji,jj)
717	zbt1 = 1.0 - zbet(ji,jj)
718	zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
719	zalf = zbt1 * zfm(ji,jj) / zpsm
720	zalf1 = 1.0 - zalf
721	ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
722	!
723	zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) )
724	zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp )
725	zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy &
726	& + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) &
727	& + ( zalf1 - zalf ) * ztemp ) )
728	zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy &
729	& + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
730	zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) )
731	zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
732	!
733	psm (ji,jj,jl) = zpsm ! optimization
734	ps0 (ji,jj,jl) = zps0
735	psx (ji,jj,jl) = zpsx
736	psxx(ji,jj,jl) = zpsxx
737	psy (ji,jj,jl) = zpsy
738	psyy(ji,jj,jl) = zpsyy
739	psxy(ji,jj,jl) = zpsxy
740	END DO
741	END DO
742
743	END DO
744	!
745	END SUBROUTINE adv_y
746
747
748	SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, &
749	& pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i )
750	!!-------------------------------------------------------------------
751	!! * ROUTINE Hbig *
752	!!
753	!! ** Purpose : Thickness correction in case advection scheme creates
754	!! abnormally tick ice or snow
755	!!
756	!! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points
757	!! (before advection) and reduce it by adapting ice concentration
758	!! 2- check whether snow thickness is larger than the surrounding 9-points
759	!! (before advection) and reduce it by sending the excess in the ocean
760	!!
761	!! ** input : Max thickness of the surrounding 9-points
762	!!-------------------------------------------------------------------
763	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
764	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max, psi_max ! max ice thick from surrounding 9-pts
765	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pes_max
766	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pei_max
767	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i
768	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
769	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i
770	!
771	INTEGER :: ji, jj, jk, jl ! dummy loop indices
772	REAL(wp) :: z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra
773	!!-------------------------------------------------------------------
774	!
775	z1_dt = 1._wp / pdt
776	!
777	DO jl = 1, jpl
778	DO jj = 1, jpj
779	DO ji = 1, jpi
780	IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
781	!
782	! ! -- check h_ip -- !
783	! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip
784	IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN
785	zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) )
786	IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN
787	pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl)
788	ENDIF
789	ENDIF
790	!
791	! ! -- check h_i -- !
792	! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i
793	zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl)
794	IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
795	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)
796	ENDIF
797	!
798	! ! -- check h_s -- !
799	! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean
800	zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl)
801	IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
802	zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )
803	!
804	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
805	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
806	!
807	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
808	pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl)
809	ENDIF
810	!
811	! ! -- check s_i -- !
812	! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean
813	zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl)
814	IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
815	zfra = psi_max(ji,jj,jl) / zsi
816	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt
817	psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra
818	ENDIF
819	!
820	ENDIF
821	END DO
822	END DO
823	END DO
824	!
825	! ! -- check e_i/v_i -- !
826	DO jl = 1, jpl
827	DO jk = 1, nlay_i
828	DO jj = 1, jpj
829	DO ji = 1, jpi
830	IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
831	! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean
832	zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl)
833	IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
834	zfra = pei_max(ji,jj,jk,jl) / zei
835	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
836	pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra
837	ENDIF
838	ENDIF
839	END DO
840	END DO
841	END DO
842	END DO
843	! ! -- check e_s/v_s -- !
844	DO jl = 1, jpl
845	DO jk = 1, nlay_s
846	DO jj = 1, jpj
847	DO ji = 1, jpi
848	IF ( pv_s(ji,jj,jl) > 0._wp ) THEN
849	! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean
850	zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl)
851	IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN
852	zfra = pes_max(ji,jj,jk,jl) / zes
853	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
854	pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra
855	ENDIF
856	ENDIF
857	END DO
858	END DO
859	END DO
860	END DO
861	!
862	END SUBROUTINE Hbig
863
864
865	SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
866	!!-------------------------------------------------------------------
867	!! * ROUTINE Hsnow *
868	!!
869	!! ** Purpose : 1- Check snow load after advection
870	!! 2- Correct pond concentration to avoid a_ip > a_i
871	!!
872	!! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface
873	!! then put the snow excess in the ocean
874	!!
875	!! ** Notes : This correction is crucial because of the call to routine icecor afterwards
876	!! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
877	!! make the snow very thick (if concentration decreases drastically)
878	!! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
879	!!-------------------------------------------------------------------
880	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
881	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip
882	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
883	!
884	INTEGER :: ji, jj, jl ! dummy loop indices
885	REAL(wp) :: z1_dt, zvs_excess, zfra
886	!!-------------------------------------------------------------------
887	!
888	z1_dt = 1._wp / pdt
889	!
890	! -- check snow load -- !
891	DO jl = 1, jpl
892	DO jj = 1, jpj
893	DO ji = 1, jpi
894	IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
895	!
896	zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
897	!
898	IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface
899	! put snow excess in the ocean
900	zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
901	wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
902	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
903	! correct snow volume and heat content
904	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
905	pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
906	ENDIF
907	!
908	ENDIF
909	END DO
910	END DO
911	END DO
912	!
913	!-- correct pond concentration to avoid a_ip > a_i -- !
914	WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
915	!
916	END SUBROUTINE Hsnow
917
918
919	SUBROUTINE adv_pra_init
920	!!-------------------------------------------------------------------
921	!! * ROUTINE adv_pra_init *
922	!!
923	!! ** Purpose : allocate and initialize arrays for Prather advection
924	!!-------------------------------------------------------------------
925	INTEGER :: ierr
926	!!-------------------------------------------------------------------
927	!
928	! !* allocate prather fields
929	ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , &
930	& sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , &
931	& sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , &
932	& sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) , &
933	& sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) , &
934	& sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , &
935	& sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , &
936	& sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) , &
937	!
938	& sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
939	& syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl) , &
940	!
941	& sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
942	& syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl) , &
943	& STAT = ierr )
944	!
945	CALL mpp_sum( 'icedyn_adv_pra', ierr )
946	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
947	!
948	CALL adv_pra_rst( 'READ' ) !* read or initialize all required files
949	!
950	END SUBROUTINE adv_pra_init
951
952
953	SUBROUTINE adv_pra_rst( cdrw, kt )
954	!!---------------------------------------------------------------------
955	!! * ROUTINE adv_pra_rst *
956	!!
957	!! ** Purpose : Read or write file in restart file
958	!!
959	!! ** Method : use of IOM library
960	!!----------------------------------------------------------------------
961	CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag
962	INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
963	!
964	INTEGER :: jk, jl ! dummy loop indices
965	INTEGER :: iter ! local integer
966	INTEGER :: id1 ! local integer
967	CHARACTER(len=25) :: znam
968	CHARACTER(len=2) :: zchar, zchar1
969	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z3d ! 3D workspace
970	!!----------------------------------------------------------------------
971	!
972	! !==========================!
973	IF( TRIM(cdrw) == 'READ' ) THEN !== Read or initialize ==!
974	! !==========================!
975	!
976	IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0
977	ELSE ; id1 = 0 ! no restart: id1=0
978	ENDIF
979	!
980	IF( id1 > 0 ) THEN ! Read the restart file !
981	!
982	! ! ice thickness
983	CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice )
984	CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice )
985	CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
986	CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
987	CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
988	! ! snow thickness
989	CALL iom_get( numrir, jpdom_autoglo, 'sxsn' , sxsn )
990	CALL iom_get( numrir, jpdom_autoglo, 'sysn' , sysn )
991	CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn )
992	CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn )
993	CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn )
994	! ! ice concentration
995	CALL iom_get( numrir, jpdom_autoglo, 'sxa' , sxa )
996	CALL iom_get( numrir, jpdom_autoglo, 'sya' , sya )
997	CALL iom_get( numrir, jpdom_autoglo, 'sxxa' , sxxa )
998	CALL iom_get( numrir, jpdom_autoglo, 'syya' , syya )
999	CALL iom_get( numrir, jpdom_autoglo, 'sxya' , sxya )
1000	! ! ice salinity
1001	CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal )
1002	CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal )
1003	CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
1004	CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
1005	CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
1006	! ! ice age
1007	CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage )
1008	CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage )
1009	CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
1010	CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
1011	CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
1012	! ! snow layers heat content
1013	DO jk = 1, nlay_s
1014	WRITE(zchar1,'(I2.2)') jk
1015	znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)
1016	znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)
1017	znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:)
1018	znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:)
1019	znam = 'sxyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:)
1020	END DO
1021	! ! ice layers heat content
1022	DO jk = 1, nlay_i
1023	WRITE(zchar1,'(I2.2)') jk
1024	znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)
1025	znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)
1026	znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:)
1027	znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:)
1028	znam = 'sxye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:)
1029	END DO
1030	!
1031	IF( ln_pnd_LEV ) THEN ! melt pond fraction
1032	IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
1033	CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap )
1034	CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap )
1035	CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
1036	CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
1037	CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
1038	! ! melt pond volume
1039	CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp )
1040	CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp )
1041	CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
1042	CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
1043	CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
1044	ELSE
1045	sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction
1046	sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume
1047	ENDIF
1048	!
1049	IF ( ln_pnd_lids ) THEN ! melt pond lid volume
1050	IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
1051	CALL iom_get( numrir, jpdom_autoglo, 'sxvl' , sxvl )
1052	CALL iom_get( numrir, jpdom_autoglo, 'syvl' , syvl )
1053	CALL iom_get( numrir, jpdom_autoglo, 'sxxvl', sxxvl )
1054	CALL iom_get( numrir, jpdom_autoglo, 'syyvl', syyvl )
1055	CALL iom_get( numrir, jpdom_autoglo, 'sxyvl', sxyvl )
1056	ELSE
1057	sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume
1058	ENDIF
1059	ENDIF
1060	ENDIF
1061	!
1062	ELSE ! start rheology from rest !
1063	!
1064	IF(lwp) WRITE(numout,*) ' ==>> start from rest OR previous run without Prather, set moments to 0'
1065	!
1066	sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness
1067	sxsn = 0._wp ; sysn = 0._wp ; sxxsn = 0._wp ; syysn = 0._wp ; sxysn = 0._wp ! snow thickness
1068	sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration
1069	sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity
1070	sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age
1071	sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content
1072	sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content
1073	IF( ln_pnd_LEV ) THEN
1074	sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction
1075	sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume
1076	IF ( ln_pnd_lids ) THEN
1077	sxvl = 0._wp; syvl = 0._wp ; sxxvl = 0._wp ; syyvl = 0._wp ; sxyvl = 0._wp ! melt pond lid volume
1078	ENDIF
1079	ENDIF
1080	ENDIF
1081	!
1082	! !=====================================!
1083	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN !== write in the ice restart file ==!
1084	! !=====================================!
1085	IF(lwp) WRITE(numout,*) '---- ice-adv-rst ----'
1086	iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1087	!
1088	!
1089	! In case Prather scheme is used for advection, write second order moments
1090	! ------------------------------------------------------------------------
1091	!
1092	! ! ice thickness
1093	CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice )
1094	CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice )
1095	CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
1096	CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
1097	CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
1098	! ! snow thickness
1099	CALL iom_rstput( iter, nitrst, numriw, 'sxsn' , sxsn )
1100	CALL iom_rstput( iter, nitrst, numriw, 'sysn' , sysn )
1101	CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn )
1102	CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn )
1103	CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn )
1104	! ! ice concentration
1105	CALL iom_rstput( iter, nitrst, numriw, 'sxa' , sxa )
1106	CALL iom_rstput( iter, nitrst, numriw, 'sya' , sya )
1107	CALL iom_rstput( iter, nitrst, numriw, 'sxxa' , sxxa )
1108	CALL iom_rstput( iter, nitrst, numriw, 'syya' , syya )
1109	CALL iom_rstput( iter, nitrst, numriw, 'sxya' , sxya )
1110	! ! ice salinity
1111	CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal )
1112	CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal )
1113	CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
1114	CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
1115	CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
1116	! ! ice age
1117	CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage )
1118	CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage )
1119	CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
1120	CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
1121	CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
1122	! ! snow layers heat content
1123	DO jk = 1, nlay_s
1124	WRITE(zchar1,'(I2.2)') jk
1125	znam = 'sxc0'//'_l'//zchar1 ; z3d(:,:,:) = sxc0 (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1126	znam = 'syc0'//'_l'//zchar1 ; z3d(:,:,:) = syc0 (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1127	znam = 'sxxc0'//'_l'//zchar1 ; z3d(:,:,:) = sxxc0(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1128	znam = 'syyc0'//'_l'//zchar1 ; z3d(:,:,:) = syyc0(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1129	znam = 'sxyc0'//'_l'//zchar1 ; z3d(:,:,:) = sxyc0(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1130	END DO
1131	! ! ice layers heat content
1132	DO jk = 1, nlay_i
1133	WRITE(zchar1,'(I2.2)') jk
1134	znam = 'sxe'//'_l'//zchar1 ; z3d(:,:,:) = sxe (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1135	znam = 'sye'//'_l'//zchar1 ; z3d(:,:,:) = sye (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1136	znam = 'sxxe'//'_l'//zchar1 ; z3d(:,:,:) = sxxe(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1137	znam = 'syye'//'_l'//zchar1 ; z3d(:,:,:) = syye(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1138	znam = 'sxye'//'_l'//zchar1 ; z3d(:,:,:) = sxye(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
1139	END DO
1140	!
1141	IF( ln_pnd_LEV ) THEN ! melt pond fraction
1142	CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap )
1143	CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap )
1144	CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
1145	CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
1146	CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
1147	! ! melt pond volume
1148	CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp )
1149	CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp )
1150	CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
1151	CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
1152	CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
1153	!
1154	IF ( ln_pnd_lids ) THEN ! melt pond lid volume
1155	CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl )
1156	CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl )
1157	CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
1158	CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
1159	CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
1160	ENDIF
1161	ENDIF
1162	!
1163	ENDIF
1164	!
1165	END SUBROUTINE adv_pra_rst
1166
1167	SUBROUTINE icemax3D( pice , pmax )
1168	!!---------------------------------------------------------------------
1169	!! * ROUTINE icemax3D *
1170	!! ** Purpose : compute the max of the 9 points around
1171	!!----------------------------------------------------------------------
1172	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input
1173	REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output
1174	REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array
1175	INTEGER :: ji, jj, jl ! dummy loop indices
1176	!!----------------------------------------------------------------------
1177	DO jl = 1, jpl
1178	DO jj = 1, jpj
1179	DO ji = 2, jpim1
1180	zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) )
1181	END DO
1182	END DO
1183	DO jj = 2, jpjm1
1184	DO ji = 2, jpim1
1185	pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
1186	END DO
1187	END DO
1188	END DO
1189	END SUBROUTINE icemax3D
1190
1191	SUBROUTINE icemax4D( pice , pmax )
1192	!!---------------------------------------------------------------------
1193	!! * ROUTINE icemax4D *
1194	!! ** Purpose : compute the max of the 9 points around
1195	!!----------------------------------------------------------------------
1196	REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input
1197	REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output
1198	REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array
1199	INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices
1200	!!----------------------------------------------------------------------
1201	jlay = SIZE( pice , 3 ) ! size of input arrays
1202	DO jl = 1, jpl
1203	DO jk = 1, jlay
1204	DO jj = 1, jpj
1205	DO ji = 2, jpim1
1206	zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) )
1207	END DO
1208	END DO
1209	DO jj = 2, jpjm1
1210	DO ji = 2, jpim1
1211	pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
1212	END DO
1213	END DO
1214	END DO
1215	END DO
1216	END SUBROUTINE icemax4D
1217
1218	#else
1219	!!----------------------------------------------------------------------
1220	!! Default option Dummy module NO SI3 sea-ice model
1221	!!----------------------------------------------------------------------
1222	#endif
1223
1224	!!======================================================================
1225	END MODULE icedyn_adv_pra

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_adv_pra.F90

Download in other formats: