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.

limrhg.F90 in branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 @ 6964

Last change on this file since 6964 was 6964, checked in by clem, 8 years ago
LIM3 rheology rewritten, see ticket #1778
Property svn:keywords set to `Id`
File size: 35.3 KB

Rev	Line
[825]	1	MODULE limrhg
	2	!!======================================================================
	3	!! * MODULE limrhg *
[834]	4	!! Ice rheology : sea ice rheology
[825]	5	!!======================================================================
[1244]	6	!! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code
	7	!! 3.0 ! 2008-03 (M. Vancoppenolle) LIM3
	8	!! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
[2528]	9	!! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas
[3680]	10	!! 3.4 ! 2011-01 (A. Porter) dynamical allocation
	11	!! 3.5 ! 2012-08 (R. Benshila) AGRIF
[6964]	12	!! 3.6 ! 2016-06 (C. Rousset) Rewriting (conserves energy)
[1244]	13	!!----------------------------------------------------------------------
[2528]	14	#if defined key_lim3 \|\| ( defined key_lim2 && ! defined key_lim2_vp )
[825]	15	!!----------------------------------------------------------------------
[2528]	16	!! 'key_lim3' OR LIM-3 sea-ice model
[2717]	17	!! 'key_lim2' AND NOT 'key_lim2_vp' EVP LIM-2 sea-ice model
[825]	18	!!----------------------------------------------------------------------
[3625]	19	!! lim_rhg : computes ice velocities
[825]	20	!!----------------------------------------------------------------------
[3625]	21	USE phycst ! Physical constant
	22	USE oce , ONLY : snwice_mass, snwice_mass_b
	23	USE par_oce ! Ocean parameters
	24	USE dom_oce ! Ocean domain
	25	USE sbc_oce ! Surface boundary condition: ocean fields
	26	USE sbc_ice ! Surface boundary condition: ice fields
[2528]	27	#if defined key_lim3
[3625]	28	USE ice ! LIM-3: ice variables
	29	USE dom_ice ! LIM-3: ice domain
	30	USE limitd_me ! LIM-3:
[2528]	31	#else
[3625]	32	USE ice_2 ! LIM-2: ice variables
	33	USE dom_ice_2 ! LIM-2: ice domain
[2528]	34	#endif
[3625]	35	USE lbclnk ! Lateral Boundary Condition / MPP link
	36	USE lib_mpp ! MPP library
	37	USE wrk_nemo ! work arrays
	38	USE in_out_manager ! I/O manager
	39	USE prtctl ! Print control
	40	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3680]	41	#if defined key_agrif && defined key_lim2
	42	USE agrif_lim2_interp
	43	#endif
[4161]	44	#if defined key_bdy
	45	USE bdyice_lim
	46	#endif
[825]	47
	48	IMPLICIT NONE
	49	PRIVATE
	50
[2715]	51	PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2)
[825]	52
[868]	53	!! * Substitutions
	54	# include "vectopt_loop_substitute.h90"
[825]	55	!!----------------------------------------------------------------------
[4161]	56	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	57	!! $Id$
[2528]	58	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	59	!!----------------------------------------------------------------------
	60	CONTAINS
	61
	62	SUBROUTINE lim_rhg( k_j1, k_jpj )
	63	!!-------------------------------------------------------------------
[834]	64	!! * SUBROUTINE lim_rhg *
	65	!! EVP-C-grid
[825]	66	!!
[834]	67	!! ** purpose : determines sea ice drift from wind stress, ice-ocean
[825]	68	!! stress and sea-surface slope. Ice-ice interaction is described by
[834]	69	!! a non-linear elasto-viscous-plastic (EVP) law including shear
	70	!! strength and a bulk rheology (Hunke and Dukowicz, 2002).
[825]	71	!!
[834]	72	!! The points in the C-grid look like this, dear reader
[825]	73	!!
[834]	74	!! (ji,jj)
	75	!! \|
	76	!! \|
	77	!! (ji-1,jj) \| (ji,jj)
	78	!! ---------
	79	!! \| \|
	80	!! \| (ji,jj) \|------(ji,jj)
	81	!! \| \|
	82	!! ---------
	83	!! (ji-1,jj-1) (ji,jj-1)
[825]	84	!!
[834]	85	!! ** Inputs : - wind forcing (stress), oceanic currents
	86	!! ice total volume (vt_i) per unit area
	87	!! snow total volume (vt_s) per unit area
[825]	88	!!
[834]	89	!! ** Action : - compute u_ice, v_ice : the components of the
	90	!! sea-ice velocity vector
	91	!! - compute delta_i, shear_i, divu_i, which are inputs
	92	!! of the ice thickness distribution
[825]	93	!!
[834]	94	!! ** Steps : 1) Compute ice snow mass, ice strength
	95	!! 2) Compute wind, oceanic stresses, mass terms and
	96	!! coriolis terms of the momentum equation
	97	!! 3) Solve the momentum equation (iterative procedure)
[6964]	98	!! 4) Recompute invariants of the strain rate tensor
[834]	99	!! which are inputs of the ITD, store stress
	100	!! for the next time step
[6964]	101	!! 5) Control prints of residual (convergence)
[834]	102	!! and charge ellipse.
	103	!! The user should make sure that the parameters
[5123]	104	!! nn_nevp, elastic time scale and rn_creepl maintain stress state
[834]	105	!! on the charge ellipse for plastic flow
	106	!! e.g. in the Canadian Archipelago
	107	!!
[6964]	108	!! ** Notes : Boundary condition for ice is chosen no-slip
	109	!! but can be adjusted with param rn_shlat
	110	!!
[2528]	111	!! References : Hunke and Dukowicz, JPO97
	112	!! Bouillon et al., Ocean Modelling 2009
	113	!!-------------------------------------------------------------------
	114	INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation
	115	INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation
[834]	116	!!
[2528]	117	INTEGER :: ji, jj ! dummy loop indices
	118	INTEGER :: jter ! local integers
[825]	119	CHARACTER (len=50) :: charout
	120
[6964]	121	REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling
	122	REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity
	123	REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013
	124	REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass
	125	REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars
	126	REAL(wp) :: zTauO, zTauE, zCor ! temporary scalars
[825]	127
[6964]	128	REAL(wp) :: zsig1, zsig2 ! internal ice stress
	129	REAL(wp) :: zresm ! Maximal error on ice velocity
	130	REAL(wp) :: zintb, zintn ! dummy argument
[3294]	131
[6964]	132	REAL(wp), POINTER, DIMENSION(:,:) :: zpresh ! temporary array for ice strength
	133	REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors
	134	REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points
	135	!
	136	REAL(wp), POINTER, DIMENSION(:,:) :: zaU , zaV ! ice fraction on U/V points
	137	REAL(wp), POINTER, DIMENSION(:,:) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points
	138	REAL(wp), POINTER, DIMENSION(:,:) :: zmf ! coriolis parameter at T points
	139	REAL(wp), POINTER, DIMENSION(:,:) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points
	140	REAL(wp), POINTER, DIMENSION(:,:) :: zspgU , zspgV ! surface pressure gradient at U/V points
	141	REAL(wp), POINTER, DIMENSION(:,:) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
	142	REAL(wp), POINTER, DIMENSION(:,:) :: zfU , zfV ! internal stresses
	143
	144	REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear
	145	REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components
	146	REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence
	147	REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope:
	148	! ocean surface (ssh_m) if ice is not embedded
	149	! ice top surface if ice is embedded
	150	REAL(wp), POINTER, DIMENSION(:,:) :: zswitchU, zswitchV ! dummy arrays
	151	REAL(wp), POINTER, DIMENSION(:,:) :: zmaskU, zmaskV ! mask for ice presence
	152	REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice
[5123]	153
[6964]	154	REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter
	155	REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity
	156	REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
[2528]	157	!!-------------------------------------------------------------------
[3294]	158
[6964]	159	CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
	160	CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
	161	CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
	162	CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
	163	CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
[3294]	164
[2528]	165	#if defined key_lim2 && ! defined key_lim2_vp
	166	# if defined key_agrif
[5123]	167	USE ice_2, vt_s => hsnm
	168	USE ice_2, vt_i => hicm
[2528]	169	# else
[5123]	170	vt_s => hsnm
	171	vt_i => hicm
[2528]	172	# endif
[5123]	173	at_i(:,:) = 1. - frld(:,:)
[2528]	174	#endif
[3680]	175	#if defined key_agrif && defined key_lim2
[5123]	176	CALL agrif_rhg_lim2_load ! First interpolation of coarse values
[3680]	177	#endif
[921]	178	!
	179	!------------------------------------------------------------------------------!
[6964]	180	! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj)
[921]	181	!------------------------------------------------------------------------------!
[6964]	182	! ocean/land mask
	183	DO jj = 1, jpjm1
	184	DO ji = 1, jpim1 ! NO vector opt.
	185	zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
	186	END DO
	187	END DO
	188	CALL lbc_lnk( zfmask, 'F', 1._wp )
[825]	189
[6964]	190	! Lateral boundary conditions on velocity (modify zfmask)
	191	zwf(:,:) = zfmask(:,:)
	192	DO jj = 2, jpjm1
	193	DO ji = fs_2, fs_jpim1 ! vector opt.
	194	IF( zfmask(ji,jj) == 0._wp ) THEN
	195	zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
	196	ENDIF
	197	END DO
	198	END DO
	199	DO jj = 2, jpjm1
	200	IF( zfmask(1,jj) == 0._wp ) THEN
	201	zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
	202	ENDIF
	203	IF( zfmask(jpi,jj) == 0._wp ) THEN
	204	zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
	205	ENDIF
	206	END DO
	207	DO ji = 2, jpim1
	208	IF( zfmask(ji,1) == 0._wp ) THEN
	209	zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
	210	ENDIF
	211	IF( zfmask(ji,jpj) == 0._wp ) THEN
	212	zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
	213	ENDIF
	214	END DO
	215	CALL lbc_lnk( zfmask, 'F', 1._wp )
	216
	217	!------------------------------------------------------------------------------!
	218	! 1) define some variables and initialize arrays
	219	!------------------------------------------------------------------------------!
	220	! ecc2: square of yield ellipse eccenticrity
	221	ecc2 = rn_ecc * rn_ecc
	222	z1_ecc2 = 1._wp / ecc2
	223
	224	! Time step for subcycling
	225	zdtevp = rdt_ice / REAL( nn_nevp )
	226	z1_dtevp = 1._wp / zdtevp
	227
	228	! alpha parameters (Bouillon 2009)
[2528]	229	#if defined key_lim3
[6964]	230	zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
	231	#else
	232	zalph1 = ( 2._wp * telast ) * z1_dtevp
[2528]	233	#endif
[6964]	234	zalph2 = zalph1 * z1_ecc2
[825]	235
[6964]	236	z1_alph1 = 1._wp / ( zalph1 + 1._wp )
	237	z1_alph2 = 1._wp / ( zalph2 + 1._wp )
	238
	239	! Initialise stress tensor
	240	zs1 (:,:) = stress1_i (:,:)
	241	zs2 (:,:) = stress2_i (:,:)
	242	zs12(:,:) = stress12_i(:,:)
	243
	244	! Ice strength
[2528]	245	#if defined key_lim3
[6964]	246	CALL lim_itd_me_icestrength( nn_icestr )
	247	zpresh(:,:) = tmask(:,:,1) * strength(:,:)
	248	#else
	249	zpresh(:,:) = tmask(:,:,1) * pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) )
[2528]	250	#endif
[825]	251
[6964]	252	! scale factors
[825]	253	DO jj = k_j1+1, k_jpj-1
[6964]	254	DO ji = fs_2, fs_jpim1
	255	z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) )
	256	z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) )
[825]	257	END DO
	258	END DO
[6964]	259
[921]	260	!
	261	!------------------------------------------------------------------------------!
	262	! 2) Wind / ocean stress, mass terms, coriolis terms
	263	!------------------------------------------------------------------------------!
	264
[3625]	265	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
[5123]	266	!
	267	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	268	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
[3625]	269	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	270	!
	271	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	272	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
[3625]	273	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
[5123]	274	!
[6964]	275	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
[5123]	276	!
[3625]	277	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	278	zpice(:,:) = ssh_m(:,:)
	279	ENDIF
	280
[825]	281	DO jj = k_j1+1, k_jpj-1
[868]	282	DO ji = fs_2, fs_jpim1
[825]	283
[6964]	284	! ice fraction at U-V points
	285	zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
	286	zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
[4990]	287
[6964]	288	! Ice/snow mass at U-V points
	289	zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )
	290	zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )
	291	zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )
	292	zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1)
	293	zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1)
[825]	294
[6964]	295	! Ocean currents at U-V points
	296	v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) &
	297	& + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
	298
	299	u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) &
	300	& + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
[825]	301
[6964]	302	! Coriolis at T points (m*f)
	303	zmf(ji,jj) = zm1 * fcor(ji,jj)
[825]	304
[6964]	305	! m/dt
	306	zmU_t(ji,jj) = zmassU * z1_dtevp
	307	zmV_t(ji,jj) = zmassV * z1_dtevp
[825]	308
[6964]	309	! Drag ice-atm.
	310	zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
	311	zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
[825]	312
[6964]	313	! Surface pressure gradient (- mgGRAD(ssh)) at U-V points
	314	zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
	315	zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
[825]	316
[6964]	317	! masks
	318	zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice
	319	zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice
[834]	320
[6964]	321	! switches
	322	zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
	323	zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
[825]	324
	325	END DO
	326	END DO
[6964]	327	CALL lbc_lnk( zmf, 'T', 1. )
[921]	328	!
	329	!------------------------------------------------------------------------------!
	330	! 3) Solution of the momentum equation, iterative procedure
	331	!------------------------------------------------------------------------------!
	332	!
[2528]	333	! !----------------------!
[5123]	334	DO jter = 1 , nn_nevp ! loop over jter !
[2528]	335	! !----------------------!
[6964]	336	IF(ln_ctl) THEN ! Convergence test
	337	DO jj = k_j1, k_jpj-1
	338	zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
	339	zv_ice(:,jj) = v_ice(:,jj)
	340	END DO
	341	ENDIF
[825]	342
[6964]	343	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
	344	DO jj = k_j1, k_jpj-1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
	345	DO ji = 1, jpim1
[825]	346
[6964]	347	! shear at F points
[5123]	348	zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
	349	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
[6964]	350	& ) * r1_e12f(ji,jj) * zfmask(ji,jj)
[825]	351
[921]	352	END DO
	353	END DO
[6964]	354	CALL lbc_lnk( zds, 'F', 1. )
[921]	355
	356	DO jj = k_j1+1, k_jpj-1
[6964]	357	DO ji = 2, jpim1 ! no vector loop
[825]	358
[6964]	359	! shear**2 at T points (doc eq. A16)
	360	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &
	361	& + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &
	362	& ) * 0.25_wp * r1_e12t(ji,jj)
	363
	364	! divergence at T points
	365	zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	366	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	367	& ) * r1_e12t(ji,jj)
	368	zdiv2 = zdiv * zdiv
	369
	370	! tension at T points
	371	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
	372	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	373	& ) * r1_e12t(ji,jj)
	374	zdt2 = zdt * zdt
	375
	376	! delta at T points
	377	zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )
[825]	378
[6964]	379	! P/delta at T points
	380	zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl )
	381
	382	! stress at T points
	383	zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
	384	zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
	385
	386	END DO
	387	END DO
	388	CALL lbc_lnk( zp_delt, 'T', 1. )
[4205]	389
[6964]	390	DO jj = k_j1, k_jpj-1
	391	DO ji = 1, jpim1
[825]	392
[6964]	393	! P/delta at F points
	394	zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
	395
	396	! stress at F points
	397	zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
[825]	398
[5123]	399	END DO
	400	END DO
[6964]	401	CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
[5429]	402
[6964]	403	! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
[921]	404	DO jj = k_j1+1, k_jpj-1
[6964]	405	DO ji = fs_2, fs_jpim1
	406
	407	! U points
	408	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
	409	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
	410	& ) * r1_e2u(ji,jj) &
	411	& + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
	412	& ) * 2._wp * r1_e1u(ji,jj) &
	413	& ) * r1_e12u(ji,jj)
	414
	415	! V points
	416	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
	417	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
	418	& ) * r1_e1v(ji,jj) &
	419	& + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
	420	& ) * 2._wp * r1_e2v(ji,jj) &
	421	& ) * r1_e12v(ji,jj)
	422
	423	! u_ice at V point
	424	u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) &
	425	& + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
	426
	427	! v_ice at U point
	428	v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) &
	429	& + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
	430
[921]	431	END DO
	432	END DO
[825]	433	!
[6964]	434	! --- Computation of ice velocity --- !
	435	! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
	436	! Bouillon et al. 2009 (eq 34-35) => stable
	437	IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
	438
[921]	439	DO jj = k_j1+1, k_jpj-1
	440	DO ji = fs_2, fs_jpim1
[825]	441
[6964]	442	! tau_io/(v_oce - v_ice)
	443	zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	444	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	445
	446	! Coriolis at V-points (energy conserving formulation)
	447	zCor = - 0.25_wp * r1_e2v(ji,jj) * &
	448	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	449	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	450
	451	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	452	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	453
	454	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	455	v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	456	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj) & ! m/dt + tau_io(only ice part)
	457	& + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	458	& ) * zmaskV(ji,jj)
[921]	459	END DO
	460	END DO
[6964]	461	CALL lbc_lnk( v_ice, 'V', -1. )
	462
[4161]	463	#if defined key_agrif && defined key_lim2
[6964]	464	CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
[3680]	465	#endif
[4333]	466	#if defined key_bdy
[6964]	467	CALL bdy_ice_lim_dyn( 'V' )
[4333]	468	#endif
[825]	469
[921]	470	DO jj = k_j1+1, k_jpj-1
	471	DO ji = fs_2, fs_jpim1
[6964]	472
	473	! tau_io/(u_oce - u_ice)
	474	zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	475	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
[834]	476
[6964]	477	! Coriolis at U-points (energy conserving formulation)
	478	zCor = 0.25_wp * r1_e1u(ji,jj) * &
	479	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	480	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	481
	482	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	483	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	484
	485	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	486	u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	487	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj) & ! m/dt + tau_io(only ice part)
	488	& + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	489	& ) * zmaskU(ji,jj)
[921]	490	END DO
	491	END DO
[6964]	492	CALL lbc_lnk( u_ice, 'U', -1. )
	493
[4161]	494	#if defined key_agrif && defined key_lim2
[6964]	495	CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
[3680]	496	#endif
[4333]	497	#if defined key_bdy
[6964]	498	CALL bdy_ice_lim_dyn( 'U' )
[4333]	499	#endif
[825]	500
[6964]	501	ELSE ! odd iterations
	502
[921]	503	DO jj = k_j1+1, k_jpj-1
	504	DO ji = fs_2, fs_jpim1
[6964]	505
	506	! tau_io/(u_oce - u_ice)
	507	zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) &
	508	& + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
[825]	509
[6964]	510	! Coriolis at U-points (energy conserving formulation)
	511	zCor = 0.25_wp * r1_e1u(ji,jj) * &
	512	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
	513	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
	514
	515	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	516	zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
	517
	518	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	519	u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	520	& ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj) & ! m/dt + tau_io(only ice part)
	521	& + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	522	& ) * zmaskU(ji,jj)
[921]	523	END DO
	524	END DO
[6964]	525	CALL lbc_lnk( u_ice, 'U', -1. )
	526
[4161]	527	#if defined key_agrif && defined key_lim2
[6964]	528	CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
[3680]	529	#endif
[4333]	530	#if defined key_bdy
[6964]	531	CALL bdy_ice_lim_dyn( 'U' )
[4333]	532	#endif
[825]	533
[6964]	534	DO jj = k_j1+1, k_jpj-1
[921]	535	DO ji = fs_2, fs_jpim1
[825]	536
[6964]	537	! tau_io/(v_oce - v_ice)
	538	zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) &
	539	& + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
	540
	541	! Coriolis at V-points (energy conserving formulation)
	542	zCor = - 0.25_wp * r1_e2v(ji,jj) * &
	543	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
	544	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
	545
	546	! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
	547	zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
	548
	549	! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
	550	v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
	551	& ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj) & ! m/dt + tau_io(only ice part)
	552	& + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin
	553	& ) * zmaskV(ji,jj)
[5123]	554	END DO
	555	END DO
[6964]	556	CALL lbc_lnk( v_ice, 'V', -1. )
	557
[4161]	558	#if defined key_agrif && defined key_lim2
[6964]	559	CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
[3680]	560	#endif
[4333]	561	#if defined key_bdy
[6964]	562	CALL bdy_ice_lim_dyn( 'V' )
[4333]	563	#endif
[825]	564
[921]	565	ENDIF
[4161]	566
[6964]	567	IF(ln_ctl) THEN ! Convergence test
	568	DO jj = k_j1+1, k_jpj-1
[5123]	569	zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
[921]	570	END DO
[5123]	571	zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
[921]	572	IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain
	573	ENDIF
[6964]	574	!
[4161]	575	! ! ==================== !
[868]	576	END DO ! end loop over jter !
[825]	577	! ! ==================== !
[921]	578	!
	579	!------------------------------------------------------------------------------!
[6964]	580	! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
[921]	581	!------------------------------------------------------------------------------!
[6964]	582	DO jj = k_j1, k_jpj-1
	583	DO ji = 1, jpim1
[866]	584
[6964]	585	! shear at F points
	586	zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) &
	587	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
	588	& ) * r1_e12f(ji,jj) * zfmask(ji,jj)
[5429]	589
[6964]	590	END DO
	591	END DO
	592	CALL lbc_lnk( zds, 'F', 1. )
	593
[868]	594	DO jj = k_j1+1, k_jpj-1
[6964]	595	DO ji = 2, jpim1 ! no vector loop
	596
	597	! tension**2 at T points
	598	zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &
	599	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
	600	& ) * r1_e12t(ji,jj)
	601	zdt2 = zdt * zdt
	602
	603	! shear**2 at T points (doc eq. A16)
	604	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &
	605	& + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &
	606	& ) * 0.25_wp * r1_e12t(ji,jj)
	607
	608	! shear at T points
	609	shear_i(ji,jj) = SQRT( zdt2 + zds2 )
[868]	610
[6964]	611	! divergence at T points
	612	divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
	613	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
	614	& ) * r1_e12t(ji,jj)
	615
	616	! delta at T points
	617	zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )
	618	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
	619	delta_i(ji,jj) = zdelta + rn_creepl * rswitch
[868]	620
[5123]	621	END DO
	622	END DO
[6964]	623	CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
	624
	625	! --- Store the stress tensor for the next time step --- !
[868]	626	stress1_i (:,:) = zs1 (:,:)
	627	stress2_i (:,:) = zs2 (:,:)
	628	stress12_i(:,:) = zs12(:,:)
	629
[921]	630	!
	631	!------------------------------------------------------------------------------!
[6964]	632	! 5) Control prints of residual and charge ellipse
[921]	633	!------------------------------------------------------------------------------!
	634	!
[834]	635	! print the residual for convergence
	636	IF(ln_ctl) THEN
[868]	637	WRITE(charout,FMT="('lim_rhg : res =',D23.16, ' iter =',I4)") zresm, jter
[834]	638	CALL prt_ctl_info(charout)
	639	CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
	640	ENDIF
[825]	641
[834]	642	! print charge ellipse
	643	! This can be desactivated once the user is sure that the stress state
	644	! lie on the charge ellipse. See Bouillon et al. 08 for more details
[825]	645	IF(ln_ctl) THEN
	646	CALL prt_ctl_info('lim_rhg : numit :',ivar1=numit)
	647	CALL prt_ctl_info('lim_rhg : nwrite :',ivar1=nwrite)
	648	CALL prt_ctl_info('lim_rhg : MOD :',ivar1=MOD(numit,nwrite))
	649	IF( MOD(numit,nwrite) .EQ. 0 ) THEN
	650	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
	651	CALL prt_ctl_info(charout)
	652	DO jj = k_j1+1, k_jpj-1
	653	DO ji = 2, jpim1
[5123]	654	IF (zpresh(ji,jj) > 1.0) THEN
[6964]	655	zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
	656	zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)*2 + 4zs12(ji,jj)2 )0.5 ) / ( 2*zpresh(ji,jj) )
[825]	657	WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
	658	CALL prt_ctl_info(charout)
	659	ENDIF
	660	END DO
	661	END DO
	662	WRITE(charout,FMT="('lim_rhg :', I4, I6, I1, I1, A10)") 2000, numit, 0, 0, ' ch. ell. '
	663	CALL prt_ctl_info(charout)
	664	ENDIF
	665	ENDIF
[2715]	666	!
[6964]	667	CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
	668	CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia )
	669	CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV )
	670	CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
	671	CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
[3294]	672
[825]	673	END SUBROUTINE lim_rhg
	674
	675	#else
	676	!!----------------------------------------------------------------------
	677	!! Default option Dummy module NO LIM sea-ice model
	678	!!----------------------------------------------------------------------
	679	CONTAINS
	680	SUBROUTINE lim_rhg( k1 , k2 ) ! Dummy routine
	681	WRITE(,) 'lim_rhg: You should not have seen this print! error?', k1, k2
	682	END SUBROUTINE lim_rhg
	683	#endif
	684
	685	!!==============================================================================
	686	END MODULE limrhg

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

Download in other formats: