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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icedyn_adv_pra.F90 in NEMO/releases/release-4.0.1/src/ICE – NEMO

Context Navigation

source: NEMO/releases/release-4.0.1/src/ICE/icedyn_adv_pra.F90 @ 11730

Last change on this file since 11730 was 11730, checked in by clem, 5 years ago
last commit: avoid division by 0 occuring on some rare occasions (CMMC and ECMWF issues). It might be machine dependant.
Property svn:keywords set to `Id`
File size: 53.8 KB

Rev	Line
[8586]	1	MODULE icedyn_adv_pra
	2	!!======================================================================
	3	!! * MODULE icedyn_adv_pra *
	4	!! sea-ice : advection => Prather scheme
	5	!!======================================================================
[9656]	6	!! History : ! 2008-03 (M. Vancoppenolle) original code
[9604]	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	8	!!--------------------------------------------------------------------
[9570]	9	#if defined key_si3
[8586]	10	!!----------------------------------------------------------------------
[9570]	11	!! 'key_si3' SI3 sea-ice model
[8586]	12	!!----------------------------------------------------------------------
[8817]	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
[8586]	17	!!----------------------------------------------------------------------
	18	USE dom_oce ! ocean domain
	19	USE ice ! sea-ice variables
	20	USE sbc_oce , ONLY : nn_fsbc ! frequency of sea-ice call
[11612]	21	USE icevar ! sea-ice: operations
[8586]	22	!
	23	USE in_out_manager ! I/O manager
	24	USE iom ! I/O manager library
	25	USE lib_mpp ! MPP library
	26	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	27	USE lbclnk ! lateral boundary conditions (or mpp links)
	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
	32	PUBLIC ice_dyn_adv_pra ! called by icedyn_adv
	33	PUBLIC adv_pra_init ! called by icedyn_adv
	34
	35	! Moments for advection
	36	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxice, syice, sxxice, syyice, sxyice ! ice thickness
	37	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsn , sysn , sxxsn , syysn , sxysn ! snow thickness
[11627]	38	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxa , sya , sxxa , syya , sxya ! ice concentration
[8586]	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxsal, sysal, sxxsal, syysal, sxysal ! ice salinity
	40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxage, syage, sxxage, syyage, sxyage ! ice age
[9271]	41	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxc0 , syc0 , sxxc0 , syyc0 , sxyc0 ! snow layers heat content
[8817]	42	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sxe , sye , sxxe , syye , sxye ! ice layers heat content
[8586]	43	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxap , syap , sxxap , syyap , sxyap ! melt pond fraction
	44	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sxvp , syvp , sxxvp , syyvp , sxyvp ! melt pond volume
[8817]	45
[8586]	46	!! * Substitutions
	47	# include "vectopt_loop_substitute.h90"
	48	!!----------------------------------------------------------------------
[9598]	49	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	50	!! $Id$
[10068]	51	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	52	!!----------------------------------------------------------------------
	53	CONTAINS
	54
	55	SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice, &
	56	& pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
	57	!!----------------------------------------------------------------------
	58	!! routine ice_dyn_adv_pra
	59	!!
	60	!! ** purpose : Computes and adds the advection trend to sea-ice
	61	!!
	62	!! ** method : Uses Prather second order scheme that advects tracers
	63	!! but also their quadratic forms. The method preserves
	64	!! tracer structures by conserving second order moments.
	65	!!
	66	!! Reference: Prather, 1986, JGR, 91, D6. 6671-6681.
	67	!!----------------------------------------------------------------------
	68	INTEGER , INTENT(in ) :: kt ! time step
	69	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity
	70	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity
	71	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
	72	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
	73	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
	74	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
	75	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
	76	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
	77	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
	78	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
	79	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
	80	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
	81	!
[11730]	82	INTEGER :: ji,jj, jk, jl, jt ! dummy loop indices
[11612]	83	INTEGER :: icycle ! number of sub-timestep for the advection
	84	REAL(wp) :: zdt ! - -
	85	REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication
[11730]	86	REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2
	87	REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx
[11612]	88	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zarea
	89	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ice, z0snw, z0ai, z0smi, z0oi
	90	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z0ap , z0vp
	91	REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) :: z0es
	92	REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) :: z0ei
[8586]	93	!!----------------------------------------------------------------------
	94	!
	95	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
	96	!
[11612]	97	! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- !
	98	! Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
	99	! this should not affect too much the stability
	100	zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
	101	zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
[8586]	102
[11612]	103	! non-blocking global communication send zcflnow and receive zcflprv
	104	CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
	105
	106	IF( zcflprv(1) > .5 ) THEN ; icycle = 2
	107	ELSE ; icycle = 1
[8586]	108	ENDIF
[11612]	109	zdt = rdt_ice / REAL(icycle)
[8586]	110
[11730]	111	! --- transport --- !
	112	zudy(:,:) = pu_ice(:,:) * e2u(:,:)
	113	zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
	114
	115	DO jt = 1, icycle
	116
	117	! record at_i before advection (for open water)
	118	zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
	119
	120	! --- transported fields --- !
	121	DO jl = 1, jpl
	122	zarea(:,:,jl) = e1e2t(:,:)
	123	z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:) ! Snow volume
	124	z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:) ! Ice volume
	125	z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:) ! Ice area
	126	z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:) ! Salt content
	127	z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:) ! Age content
	128	DO jk = 1, nlay_s
	129	z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
	130	END DO
	131	DO jk = 1, nlay_i
	132	z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content
	133	END DO
	134	IF ( ln_pnd_H12 ) THEN
	135	z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction
	136	z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume
	137	ENDIF
[9271]	138	END DO
[11730]	139	!
	140	! !--------------------------------------------!
	141	IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==!
	142	! !--------------------------------------------!
	143	CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
	144	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
	145	CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
	146	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
	147	CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
	148	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
	149	CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
	150	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
	151	CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
	152	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
[11612]	153	!
[11730]	154	DO jk = 1, nlay_s !--- snow heat content
	155	CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
	156	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
	157	CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
	158	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
[8586]	159	END DO
[11730]	160	DO jk = 1, nlay_i !--- ice heat content
	161	CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
	162	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
	163	CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
	164	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
[11612]	165	END DO
	166	!
	167	IF ( ln_pnd_H12 ) THEN
[11730]	168	CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
	169	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
	170	CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
	171	CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
[11612]	172	ENDIF
[11730]	173	! !--------------------------------------------!
	174	ELSE !== even ice time step: adv_y then adv_x ==!
	175	! !--------------------------------------------!
	176	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
	177	CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
	178	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn ) !--- snow volume
	179	CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn , sxxsn , sysn , syysn , sxysn )
	180	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
	181	CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
	182	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya ) !--- ice concentration
	183	CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai , sxa , sxxa , sya , syya , sxya )
	184	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
	185	CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi , sxage , sxxage , syage , syyage , sxyage )
	186	DO jk = 1, nlay_s !--- snow heat content
	187	CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
	188	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
	189	CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:), &
	190	& sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
[8586]	191	END DO
[11730]	192	DO jk = 1, nlay_i !--- ice heat content
	193	CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
	194	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
	195	CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:), &
	196	& sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
[11612]	197	END DO
	198	IF ( ln_pnd_H12 ) THEN
[11730]	199	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction
	200	CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
	201	CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) !--- melt pond volume
	202	CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
[11612]	203	ENDIF
[11730]	204	!
	205	ENDIF
[8586]	206
[11730]	207	! --- Recover the properties from their contents --- !
	208	DO jl = 1, jpl
	209	pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	210	pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	211	psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	212	poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	213	pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	214	DO jk = 1, nlay_s
	215	pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	216	END DO
	217	DO jk = 1, nlay_i
	218	pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	219	END DO
	220	IF ( ln_pnd_H12 ) THEN
	221	pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	222	pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
	223	ENDIF
[9271]	224	END DO
[11730]	225	!
	226	! derive open water from ice concentration
	227	zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
	228	DO jj = 2, jpjm1
	229	DO ji = fs_2, fs_jpim1
	230	pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water
	231	& - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
	232	END DO
[8586]	233	END DO
[11730]	234	CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. )
	235	!
	236	! --- Ensure non-negative fields --- !
	237	! Remove negative values (conservation is ensured)
	238	! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
	239	CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
	240	!
	241	! --- Ensure snow load is not too big --- !
	242	CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
	243	!
[8586]	244	END DO
	245	!
	246	IF( lrst_ice ) CALL adv_pra_rst( 'WRITE', kt ) !* write Prather fields in the restart file
	247	!
	248	END SUBROUTINE ice_dyn_adv_pra
	249
[8817]	250
[11612]	251	SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 , &
[8586]	252	& psx, psxx, psy , psyy, psxy )
	253	!!----------------------------------------------------------------------
	254	!! routine adv_x
	255	!!
	256	!! ** purpose : Computes and adds the advection trend to sea-ice
	257	!! variable on x axis
	258	!!----------------------------------------------------------------------
[11612]	259	REAL(wp) , INTENT(in ) :: pdt ! the time step
	260	REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
	261	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: put ! i-direction ice velocity at U-point [m/s]
	262	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
	263	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
	264	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
	265	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
[8586]	266	!!
[11612]	267	INTEGER :: ji, jj, jl, jcat ! dummy loop indices
	268	REAL(wp) :: zs1max, zslpmax, ztemp ! local scalars
[8586]	269	REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
	270	REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
	271	REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace
	272	REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - -
	273	REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
	274	!-----------------------------------------------------------------------
[11612]	275	!
	276	jcat = SIZE( ps0 , 3 ) ! size of input arrays
	277	!
	278	DO jl = 1, jcat ! loop on categories
	279	!
	280	! Limitation of moments.
	281	DO jj = 2, jpjm1
	282	DO ji = 1, jpi
	283	! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
	284	psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
	285	!
	286	zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
	287	zs1max = 1.5 * zslpmax
	288	zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
	289	zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), &
	290	& MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) )
	291	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
[8586]	292
[11612]	293	ps0 (ji,jj,jl) = zslpmax
	294	psx (ji,jj,jl) = zs1new * rswitch
	295	psxx(ji,jj,jl) = zs2new * rswitch
	296	psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
	297	psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
	298	psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
	299	END DO
	300	END DO
[8586]	301
[11612]	302	! Calculate fluxes and moments between boxes i<-->i+1
	303	DO jj = 2, jpjm1 ! Flux from i to i+1 WHEN u GT 0
	304	DO ji = 1, jpi
	305	zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
[11730]	306	zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
[11612]	307	zalfq = zalf * zalf
	308	zalf1 = 1.0 - zalf
	309	zalf1q = zalf1 * zalf1
	310	!
	311	zfm (ji,jj) = zalf * psm (ji,jj,jl)
	312	zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
	313	zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
	314	zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq
	315	zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
	316	zfxy(ji,jj) = zalfq * psxy(ji,jj,jl)
	317	zfyy(ji,jj) = zalf * psyy(ji,jj,jl)
[8586]	318
[11612]	319	! Readjust moments remaining in the box.
	320	psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj)
	321	ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj)
	322	psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
	323	psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl)
	324	psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj)
	325	psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj)
	326	psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl)
	327	END DO
[8586]	328	END DO
	329
[11612]	330	DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0.
	331	DO ji = 1, fs_jpim1
[11730]	332	zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)
[11612]	333	zalg (ji,jj) = zalf
	334	zalfq = zalf * zalf
	335	zalf1 = 1.0 - zalf
	336	zalg1 (ji,jj) = zalf1
	337	zalf1q = zalf1 * zalf1
	338	zalg1q(ji,jj) = zalf1q
	339	!
	340	zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl)
	341	zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) &
	342	& - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
	343	zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
	344	zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq
	345	zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
	346	zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl)
	347	zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl)
	348	END DO
[8586]	349	END DO
	350
[11612]	351	DO jj = 2, jpjm1 ! Readjust moments remaining in the box.
	352	DO ji = fs_2, fs_jpim1
	353	zbt = zbet(ji-1,jj)
	354	zbt1 = 1.0 - zbet(ji-1,jj)
	355	!
	356	psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
	357	ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
	358	psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
	359	psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
	360	psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
	361	psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
	362	psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
	363	END DO
[8586]	364	END DO
	365
[11612]	366	! Put the temporary moments into appropriate neighboring boxes.
	367	DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0.
	368	DO ji = fs_2, fs_jpim1
	369	zbt = zbet(ji-1,jj)
	370	zbt1 = 1.0 - zbet(ji-1,jj)
	371	psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
	372	zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
	373	zalf1 = 1.0 - zalf
	374	ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
	375	!
	376	ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
	377	psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
	378	psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
	379	& + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
	380	& + zbt1 * psxx(ji,jj,jl)
	381	psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) &
	382	& + 3.0 * (- zalf1zfy(ji-1,jj) + zalf psy(ji,jj,jl) ) ) &
	383	& + zbt1 * psxy(ji,jj,jl)
	384	psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
	385	psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
	386	END DO
[8586]	387	END DO
	388
[11612]	389	DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0.
	390	DO ji = fs_2, fs_jpim1
	391	zbt = zbet(ji,jj)
	392	zbt1 = 1.0 - zbet(ji,jj)
	393	psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
	394	zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
	395	zalf1 = 1.0 - zalf
	396	ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
	397	!
	398	ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
	399	psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
	400	psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
	401	& + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) &
	402	& + ( zalf1 - zalf ) * ztemp ) )
	403	psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) &
	404	& + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
	405	psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
	406	psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
	407	END DO
[8586]	408	END DO
	409
	410	END DO
	411
	412	!-- Lateral boundary conditions
[11612]	413	CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. &
	414	& , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes
	415	& , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. )
[8586]	416	!
	417	END SUBROUTINE adv_x
	418
	419
[11612]	420	SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 , &
[8586]	421	& psx, psxx, psy , psyy, psxy )
	422	!!---------------------------------------------------------------------
	423	!! routine adv_y
	424	!!
	425	!! ** purpose : Computes and adds the advection trend to sea-ice
	426	!! variable on y axis
	427	!!---------------------------------------------------------------------
[11612]	428	REAL(wp) , INTENT(in ) :: pdt ! time step
	429	REAL(wp) , INTENT(in ) :: pcrh ! call adv_x then adv_y (=1) or the opposite (=0)
	430	REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvt ! j-direction ice velocity at V-point [m/s]
	431	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psm ! area
	432	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ps0 ! field to be advected
	433	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psx , psy ! 1st moments
	434	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: psxx, psyy, psxy ! 2nd moments
[8586]	435	!!
[11612]	436	INTEGER :: ji, jj, jl, jcat ! dummy loop indices
	437	REAL(wp) :: zs1max, zslpmax, ztemp ! temporary scalars
[8586]	438	REAL(wp) :: zs1new, zalf , zalfq , zbt ! - -
	439	REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - -
	440	REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace
	441	REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - -
	442	REAL(wp), DIMENSION(jpi,jpj) :: zalg, zalg1, zalg1q ! - -
	443	!---------------------------------------------------------------------
[11612]	444	!
	445	jcat = SIZE( ps0 , 3 ) ! size of input arrays
	446	!
	447	DO jl = 1, jcat ! loop on categories
	448	!
	449	! Limitation of moments.
	450	DO jj = 1, jpj
	451	DO ji = fs_2, fs_jpim1
	452	! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
	453	psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
	454	!
	455	zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
	456	zs1max = 1.5 * zslpmax
	457	zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
	458	zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), &
	459	& MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) )
	460	rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask
	461	!
	462	ps0 (ji,jj,jl) = zslpmax
	463	psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
	464	psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
	465	psy (ji,jj,jl) = zs1new * rswitch
	466	psyy(ji,jj,jl) = zs2new * rswitch
	467	psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
	468	END DO
[8586]	469	END DO
[11612]	470
	471	! Calculate fluxes and moments between boxes j<-->j+1
	472	DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0
	473	DO ji = fs_2, fs_jpim1
	474	zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
[11730]	475	zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
[11612]	476	zalfq = zalf * zalf
	477	zalf1 = 1.0 - zalf
	478	zalf1q = zalf1 * zalf1
	479	!
	480	zfm (ji,jj) = zalf * psm(ji,jj,jl)
	481	zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) )
	482	zfy (ji,jj) = zalfq ( psy(ji,jj,jl) + 3.0zalf1*psyy(ji,jj,jl) )
	483	zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl)
	484	zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
	485	zfxy(ji,jj) = zalfq * psxy(ji,jj,jl)
	486	zfxx(ji,jj) = zalf * psxx(ji,jj,jl)
	487	!
	488	! Readjust moments remaining in the box.
	489	psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj)
	490	ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj)
	491	psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
	492	psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl)
	493	psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj)
	494	psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj)
	495	psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl)
	496	END DO
[8586]	497	END DO
[11612]	498	!
	499	DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0.
	500	DO ji = fs_2, fs_jpim1
[11730]	501	zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)
[11612]	502	zalg (ji,jj) = zalf
	503	zalfq = zalf * zalf
	504	zalf1 = 1.0 - zalf
	505	zalg1 (ji,jj) = zalf1
	506	zalf1q = zalf1 * zalf1
	507	zalg1q(ji,jj) = zalf1q
	508	!
	509	zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl)
	510	zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) &
	511	& - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
	512	zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
	513	zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq
	514	zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
	515	zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl)
	516	zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl)
	517	END DO
[8586]	518	END DO
	519
[11612]	520	! Readjust moments remaining in the box.
	521	DO jj = 2, jpjm1
	522	DO ji = fs_2, fs_jpim1
	523	zbt = zbet(ji,jj-1)
	524	zbt1 = ( 1.0 - zbet(ji,jj-1) )
	525	!
	526	psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
	527	ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
	528	psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
	529	psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
	530	psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
	531	psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
	532	psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
	533	END DO
[8586]	534	END DO
	535
[11612]	536	! Put the temporary moments into appropriate neighboring boxes.
	537	DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0.
	538	DO ji = fs_2, fs_jpim1
	539	zbt = zbet(ji,jj-1)
	540	zbt1 = 1.0 - zbet(ji,jj-1)
	541	psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)
	542	zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)
	543	zalf1 = 1.0 - zalf
	544	ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
	545	!
	546	ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
	547	psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) &
	548	& + zbt1 * psy(ji,jj,jl)
	549	psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) &
	550	& + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &
	551	& + zbt1 * psyy(ji,jj,jl)
	552	psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) &
	553	& + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) &
	554	& + zbt1 * psxy(ji,jj,jl)
	555	psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
	556	psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
	557	END DO
[8586]	558	END DO
	559
[11612]	560	DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0.
	561	DO ji = fs_2, fs_jpim1
	562	zbt = zbet(ji,jj)
	563	zbt1 = 1.0 - zbet(ji,jj)
	564	psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
	565	zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
	566	zalf1 = 1.0 - zalf
	567	ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
	568	!
	569	ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
	570	psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
	571	psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
	572	& + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) &
	573	& + ( zalf1 - zalf ) * ztemp ) )
	574	psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) &
	575	& + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
	576	psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
	577	psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
	578	END DO
[8586]	579	END DO
[11612]	580
[8586]	581	END DO
	582
	583	!-- Lateral boundary conditions
[11612]	584	CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T', 1., ps0 , 'T', 1. &
	585	& , psx , 'T', -1., psy , 'T', -1. & ! caution gradient ==> the sign changes
	586	& , psxx , 'T', 1., psyy, 'T', 1. , psxy, 'T', 1. )
[8586]	587	!
	588	END SUBROUTINE adv_y
	589
[8817]	590
[11632]	591	SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
	592	!!-------------------------------------------------------------------
	593	!! * ROUTINE Hsnow *
	594	!!
	595	!! ** Purpose : 1- Check snow load after advection
	596	!! 2- Correct pond concentration to avoid a_ip > a_i
	597	!!
	598	!! ** Method : If snow load makes snow-ice interface to deplet below the ocean surface
	599	!! then put the snow excess in the ocean
	600	!!
	601	!! ** Notes : This correction is crucial because of the call to routine icecor afterwards
	602	!! which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
	603	!! make the snow very thick (if concentration decreases drastically)
	604	!! This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
	605	!!-------------------------------------------------------------------
	606	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
	607	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, pa_i, pa_ip
	608	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s
	609	!
	610	INTEGER :: ji, jj, jl ! dummy loop indices
	611	REAL(wp) :: z1_dt, zvs_excess, zfra
	612	!!-------------------------------------------------------------------
	613	!
	614	z1_dt = 1._wp / pdt
	615	!
	616	! -- check snow load -- !
	617	DO jl = 1, jpl
	618	DO jj = 1, jpj
	619	DO ji = 1, jpi
	620	IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
	621	!
	622	zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
	623	!
	624	IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface
	625	! put snow excess in the ocean
	626	zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
	627	wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
	628	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
	629	! correct snow volume and heat content
	630	pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
	631	pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess
	632	ENDIF
	633	!
	634	ENDIF
	635	END DO
	636	END DO
	637	END DO
	638	!
	639	!-- correct pond concentration to avoid a_ip > a_i -- !
	640	WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:)
	641	!
	642	END SUBROUTINE Hsnow
	643
	644
[8586]	645	SUBROUTINE adv_pra_init
	646	!!-------------------------------------------------------------------
	647	!! * ROUTINE adv_pra_init *
	648	!!
	649	!! ** Purpose : allocate and initialize arrays for Prather advection
	650	!!-------------------------------------------------------------------
	651	INTEGER :: ierr
	652	!!-------------------------------------------------------------------
[8817]	653	!
	654	! !* allocate prather fields
[11730]	655	ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) , &
[8586]	656	& sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) , &
	657	& sxa (jpi,jpj,jpl) , sya (jpi,jpj,jpl) , sxxa (jpi,jpj,jpl) , syya (jpi,jpj,jpl) , sxya (jpi,jpj,jpl) , &
	658	& sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) , &
	659	& sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) , &
	660	& sxap(jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) , &
	661	& sxvp(jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) , &
[9271]	662	!
	663	& sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
	664	& syyc0(jpi,jpj,nlay_s,jpl) , sxyc0(jpi,jpj,nlay_s,jpl) , &
	665	!
	666	& sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe (jpi,jpj,nlay_i,jpl) , &
	667	& syye (jpi,jpj,nlay_i,jpl) , sxye (jpi,jpj,nlay_i,jpl) , &
[8586]	668	& STAT = ierr )
	669	!
[10425]	670	CALL mpp_sum( 'icedyn_adv_pra', ierr )
[8586]	671	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'adv_pra_init : unable to allocate ice arrays for Prather advection scheme')
	672	!
[8817]	673	CALL adv_pra_rst( 'READ' ) !* read or initialize all required files
[8586]	674	!
	675	END SUBROUTINE adv_pra_init
	676
[8817]	677
[8586]	678	SUBROUTINE adv_pra_rst( cdrw, kt )
	679	!!---------------------------------------------------------------------
	680	!! * ROUTINE adv_pra_rst *
	681	!!
[11612]	682	!! ** Purpose : Read or write file in restart file
[8586]	683	!!
	684	!! ** Method : use of IOM library
	685	!!----------------------------------------------------------------------
	686	CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	687	INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step
	688	!
	689	INTEGER :: jk, jl ! dummy loop indices
	690	INTEGER :: iter ! local integer
	691	INTEGER :: id1 ! local integer
	692	CHARACTER(len=25) :: znam
	693	CHARACTER(len=2) :: zchar, zchar1
[8817]	694	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z3d ! 3D workspace
[8586]	695	!!----------------------------------------------------------------------
	696	!
[8817]	697	! !==========================!
	698	IF( TRIM(cdrw) == 'READ' ) THEN !== Read or initialize ==!
	699	! !==========================!
	700	!
[11730]	701	IF( ln_rstart ) THEN ; id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. ) ! file exist: id1>0
[8817]	702	ELSE ; id1 = 0 ! no restart: id1=0
	703	ENDIF
	704	!
	705	IF( id1 > 0 ) THEN ! Read the restart file !
[8586]	706	!
[8817]	707	! ! ice thickness
	708	CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice )
	709	CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice )
	710	CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
	711	CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
	712	CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
	713	! ! snow thickness
	714	CALL iom_get( numrir, jpdom_autoglo, 'sxsn' , sxsn )
	715	CALL iom_get( numrir, jpdom_autoglo, 'sysn' , sysn )
	716	CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn )
	717	CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn )
	718	CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn )
[11627]	719	! ! ice concentration
[8817]	720	CALL iom_get( numrir, jpdom_autoglo, 'sxa' , sxa )
	721	CALL iom_get( numrir, jpdom_autoglo, 'sya' , sya )
	722	CALL iom_get( numrir, jpdom_autoglo, 'sxxa' , sxxa )
	723	CALL iom_get( numrir, jpdom_autoglo, 'syya' , syya )
	724	CALL iom_get( numrir, jpdom_autoglo, 'sxya' , sxya )
	725	! ! ice salinity
	726	CALL iom_get( numrir, jpdom_autoglo, 'sxsal' , sxsal )
	727	CALL iom_get( numrir, jpdom_autoglo, 'sysal' , sysal )
	728	CALL iom_get( numrir, jpdom_autoglo, 'sxxsal', sxxsal )
	729	CALL iom_get( numrir, jpdom_autoglo, 'syysal', syysal )
	730	CALL iom_get( numrir, jpdom_autoglo, 'sxysal', sxysal )
	731	! ! ice age
	732	CALL iom_get( numrir, jpdom_autoglo, 'sxage' , sxage )
	733	CALL iom_get( numrir, jpdom_autoglo, 'syage' , syage )
	734	CALL iom_get( numrir, jpdom_autoglo, 'sxxage', sxxage )
	735	CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
	736	CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
[9271]	737	! ! snow layers heat content
	738	DO jk = 1, nlay_s
	739	WRITE(zchar1,'(I2.2)') jk
	740	znam = 'sxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxc0 (:,:,jk,:) = z3d(:,:,:)
	741	znam = 'syc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; syc0 (:,:,jk,:) = z3d(:,:,:)
	742	znam = 'sxxc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxxc0(:,:,jk,:) = z3d(:,:,:)
	743	znam = 'syyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; syyc0(:,:,jk,:) = z3d(:,:,:)
	744	znam = 'sxyc0'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxyc0(:,:,jk,:) = z3d(:,:,:)
	745	END DO
[8817]	746	! ! ice layers heat content
[9271]	747	DO jk = 1, nlay_i
[8817]	748	WRITE(zchar1,'(I2.2)') jk
[9123]	749	znam = 'sxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxe (:,:,jk,:) = z3d(:,:,:)
	750	znam = 'sye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sye (:,:,jk,:) = z3d(:,:,:)
	751	znam = 'sxxe'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxxe(:,:,jk,:) = z3d(:,:,:)
	752	znam = 'syye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; syye(:,:,jk,:) = z3d(:,:,:)
	753	znam = 'sxye'//'_l'//zchar1 ; CALL iom_get( numrir, jpdom_autoglo, znam , z3d ) ; sxye(:,:,jk,:) = z3d(:,:,:)
[8817]	754	END DO
[8586]	755	!
[8817]	756	IF( ln_pnd_H12 ) THEN ! melt pond fraction
	757	CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap )
	758	CALL iom_get( numrir, jpdom_autoglo, 'syap' , syap )
	759	CALL iom_get( numrir, jpdom_autoglo, 'sxxap', sxxap )
	760	CALL iom_get( numrir, jpdom_autoglo, 'syyap', syyap )
	761	CALL iom_get( numrir, jpdom_autoglo, 'sxyap', sxyap )
	762	! ! melt pond volume
	763	CALL iom_get( numrir, jpdom_autoglo, 'sxvp' , sxvp )
	764	CALL iom_get( numrir, jpdom_autoglo, 'syvp' , syvp )
	765	CALL iom_get( numrir, jpdom_autoglo, 'sxxvp', sxxvp )
	766	CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
	767	CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
[8586]	768	ENDIF
	769	!
[8817]	770	ELSE ! start rheology from rest !
	771	!
	772	IF(lwp) WRITE(numout,*) ' ==>> start from rest OR previous run without Prather, set moments to 0'
	773	!
	774	sxice = 0._wp ; syice = 0._wp ; sxxice = 0._wp ; syyice = 0._wp ; sxyice = 0._wp ! ice thickness
	775	sxsn = 0._wp ; sysn = 0._wp ; sxxsn = 0._wp ; syysn = 0._wp ; sxysn = 0._wp ! snow thickness
[11627]	776	sxa = 0._wp ; sya = 0._wp ; sxxa = 0._wp ; syya = 0._wp ; sxya = 0._wp ! ice concentration
[8817]	777	sxsal = 0._wp ; sysal = 0._wp ; sxxsal = 0._wp ; syysal = 0._wp ; sxysal = 0._wp ! ice salinity
	778	sxage = 0._wp ; syage = 0._wp ; sxxage = 0._wp ; syyage = 0._wp ; sxyage = 0._wp ! ice age
[9271]	779	sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content
[8817]	780	sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content
	781	IF( ln_pnd_H12 ) THEN
	782	sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction
	783	sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume
[8586]	784	ENDIF
	785	ENDIF
	786	!
[8817]	787	! !=====================================!
	788	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN !== write in the ice restart file ==!
	789	! !=====================================!
	790	IF(lwp) WRITE(numout,*) '---- ice-adv-rst ----'
[8586]	791	iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1
	792	!
[8817]	793	!
	794	! In case Prather scheme is used for advection, write second order moments
	795	! ------------------------------------------------------------------------
	796	!
	797	! ! ice thickness
	798	CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice )
	799	CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice )
	800	CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
	801	CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
	802	CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
	803	! ! snow thickness
	804	CALL iom_rstput( iter, nitrst, numriw, 'sxsn' , sxsn )
	805	CALL iom_rstput( iter, nitrst, numriw, 'sysn' , sysn )
	806	CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn )
	807	CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn )
	808	CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn )
[11627]	809	! ! ice concentration
[8817]	810	CALL iom_rstput( iter, nitrst, numriw, 'sxa' , sxa )
	811	CALL iom_rstput( iter, nitrst, numriw, 'sya' , sya )
	812	CALL iom_rstput( iter, nitrst, numriw, 'sxxa' , sxxa )
	813	CALL iom_rstput( iter, nitrst, numriw, 'syya' , syya )
	814	CALL iom_rstput( iter, nitrst, numriw, 'sxya' , sxya )
	815	! ! ice salinity
	816	CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal )
	817	CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal )
	818	CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
	819	CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
	820	CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
	821	! ! ice age
	822	CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage )
	823	CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage )
	824	CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
	825	CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
	826	CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
[9271]	827	! ! snow layers heat content
	828	DO jk = 1, nlay_s
	829	WRITE(zchar1,'(I2.2)') jk
	830	znam = 'sxc0'//'_l'//zchar1 ; z3d(:,:,:) = sxc0 (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	831	znam = 'syc0'//'_l'//zchar1 ; z3d(:,:,:) = syc0 (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	832	znam = 'sxxc0'//'_l'//zchar1 ; z3d(:,:,:) = sxxc0(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	833	znam = 'syyc0'//'_l'//zchar1 ; z3d(:,:,:) = syyc0(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	834	znam = 'sxyc0'//'_l'//zchar1 ; z3d(:,:,:) = sxyc0(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	835	END DO
[8817]	836	! ! ice layers heat content
[9271]	837	DO jk = 1, nlay_i
[8817]	838	WRITE(zchar1,'(I2.2)') jk
[9123]	839	znam = 'sxe'//'_l'//zchar1 ; z3d(:,:,:) = sxe (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	840	znam = 'sye'//'_l'//zchar1 ; z3d(:,:,:) = sye (:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	841	znam = 'sxxe'//'_l'//zchar1 ; z3d(:,:,:) = sxxe(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	842	znam = 'syye'//'_l'//zchar1 ; z3d(:,:,:) = syye(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
	843	znam = 'sxye'//'_l'//zchar1 ; z3d(:,:,:) = sxye(:,:,jk,:) ; CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
[8586]	844	END DO
[8817]	845	!
	846	IF( ln_pnd_H12 ) THEN ! melt pond fraction
	847	CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap )
	848	CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap )
	849	CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
	850	CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
	851	CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
	852	! ! melt pond volume
	853	CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp )
	854	CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp )
	855	CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
	856	CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
	857	CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
[8586]	858	ENDIF
	859	!
	860	ENDIF
	861	!
	862	END SUBROUTINE adv_pra_rst
	863
	864	#else
	865	!!----------------------------------------------------------------------
[9570]	866	!! Default option Dummy module NO SI3 sea-ice model
[8586]	867	!!----------------------------------------------------------------------
	868	#endif
	869
	870	!!======================================================================
	871	END MODULE icedyn_adv_pra

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

Download in other formats: