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.

icerdgrft.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90 @ 8517

Last change on this file since 8517 was 8517, checked in by clem, 7 years ago
changes in style - part6 - one more round
File size: 49.3 KB

Rev	Line
[8409]	1	MODULE icerdgrft
	2	!!======================================================================
	3	!! * MODULE icerdgrft *
	4	!! LIM-3 : Mechanical impact on ice thickness distribution
	5	!!======================================================================
	6	!! History : LIM ! 2006-02 (M. Vancoppenolle) Original code
	7	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
	8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	9	!!----------------------------------------------------------------------
	10	#if defined key_lim3
	11	!!----------------------------------------------------------------------
	12	!! 'key_lim3' LIM-3 sea-ice model
	13	!!----------------------------------------------------------------------
[8486]	14	USE par_oce ! ocean parameters
	15	USE dom_oce ! ocean domain
	16	USE phycst ! physical constants (ocean directory)
	17	USE sbc_oce , ONLY : sss_m, sst_m ! surface boundary condition: ocean fields
	18	USE ice1D ! sea-ice: thermodynamics
	19	USE ice ! sea-ice: variables
	20	USE icevar ! sea-ice: operations
	21	USE icectl ! sea-ice: control prints
[8409]	22	!
[8486]	23	USE lbclnk ! lateral boundary condition - MPP exchanges
	24	USE lib_mpp ! MPP library
	25	USE in_out_manager ! I/O manager
	26	USE iom ! I/O manager
	27	USE lib_fortran ! glob_sum
	28	USE timing ! Timing
[8409]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
	33	PUBLIC ice_rdgrft ! called by ice_stp
[8512]	34	PUBLIC ice_rdgrft_strength ! called by icerhg_evp
[8486]	35	PUBLIC ice_rdgrft_init ! called by ice_stp
	36	PUBLIC ice_rdgrft_alloc ! called by ice_init
[8409]	37
	38	! Variables shared among ridging subroutines
	39	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area
	40	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged
	41	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/closing associated w/ category n
	42	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness
	43	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness
	44	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice
	45	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! thickness of ridging ice / mean ridge thickness
	46	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging
	47	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting
[8498]	48	!
[8409]	49	REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier
	50	REAL(wp), PARAMETER :: kraft = 0.5_wp ! rafting multipliyer
[8498]	51	REAL(wp) :: zdrho !
[8409]	52	!
[8512]	53	! namelist (namice_rdgrft)
[8515]	54	REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging
[8512]	55	LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975))
	56	REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging
	57	LOGICAL :: ln_partf_exp ! participation function exponential (Lipscomb et al. (2007))
	58	REAL(wp) :: rn_astar ! equivalent of G* for an exponential participation function
	59	LOGICAL :: ln_ridging ! ridging of ice or not
	60	REAL(wp) :: rn_hstar ! thickness that determines the maximal thickness of ridged ice
	61	REAL(wp) :: rn_porordg ! initial porosity of ridges (0.3 regular value)
	62	REAL(wp) :: rn_fsnwrdg ! fractional snow loss to the ocean during ridging
	63	REAL(wp) :: rn_fpndrdg ! fractional pond loss to the ocean during ridging
	64	LOGICAL :: ln_rafting ! rafting of ice or not
	65	REAL(wp) :: rn_hraft ! threshold thickness (m) for rafting / ridging
	66	REAL(wp) :: rn_craft ! coefficient for smoothness of the hyperbolic tangent in rafting
	67	REAL(wp) :: rn_fsnwrft ! fractional snow loss to the ocean during rafting
	68	REAL(wp) :: rn_fpndrft ! fractional pond loss to the ocean during rafting
[8498]	69	!
[8409]	70	!!----------------------------------------------------------------------
[8486]	71	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8409]	72	!! $Id: icerdgrft.F90 8378 2017-07-26 13:55:59Z clem $
	73	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	74	!!----------------------------------------------------------------------
	75	CONTAINS
	76
	77	INTEGER FUNCTION ice_rdgrft_alloc()
	78	!!---------------------------------------------------------------------!
	79	!! * ROUTINE ice_rdgrft_alloc *
	80	!!---------------------------------------------------------------------!
[8486]	81	ALLOCATE( asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj) , &
	82	& hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , &
	83	& hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
[8409]	84	!
	85	IF( ice_rdgrft_alloc /= 0 ) CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' )
	86	!
	87	END FUNCTION ice_rdgrft_alloc
	88
	89
[8426]	90	SUBROUTINE ice_rdgrft( kt )
[8409]	91	!!---------------------------------------------------------------------!
	92	!! * ROUTINE ice_rdgrft *
	93	!!
	94	!! ** Purpose : computes the mechanical redistribution of ice thickness
	95	!!
	96	!! ** Method : Steps :
	97	!! 1) Thickness categories boundaries, ice / o.w. concentrations
	98	!! Ridge preparation
	99	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
	100	!! 3) Ridging iteration
	101	!! 4) Ridging diagnostics
	102	!! 5) Heat, salt and freshwater fluxes
	103	!! 6) Compute increments of tate variables and come back to old values
	104	!!
	105	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
	106	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
	107	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
	108	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
	109	!! Bitz et al., JGR, 2001
	110	!! Amundrud and Melling, JGR 2005
	111	!! Babko et al., JGR 2002
	112	!!
	113	!! This routine is based on CICE code and authors William H. Lipscomb,
	114	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
	115	!!--------------------------------------------------------------------!
[8426]	116	INTEGER, INTENT(in) :: kt ! number of iteration
	117	!!
[8486]	118	INTEGER :: ji, jj, jk, jl ! dummy loop index
	119	INTEGER :: niter ! local integer
	120	INTEGER :: iterate_ridging ! if =1, repeat the ridging
[8498]	121	REAL(wp) :: za ! local scalar
[8486]	122	REAL(wp), DIMENSION(jpi,jpj) :: closing_net ! net rate at which area is removed (1/s)
	123	! ! (ridging ice area - area of new ridges) / dt
	124	REAL(wp), DIMENSION(jpi,jpj) :: divu_adv ! divu as implied by transport scheme (1/s)
	125	REAL(wp), DIMENSION(jpi,jpj) :: opning ! rate of opening due to divergence/shear
	126	REAL(wp), DIMENSION(jpi,jpj) :: closing_gross ! rate at which area removed, not counting area of new ridges
[8409]	127	!
	128	INTEGER, PARAMETER :: nitermax = 20
	129	!!-----------------------------------------------------------------------------
[8517]	130	! controls
	131	IF( nn_timing == 1 ) CALL timing_start('icerdgrft') ! timing
	132	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[8409]	133
[8498]	134	IF( kt == nit000 ) THEN
	135	IF(lwp) WRITE(numout,*)
[8512]	136	IF(lwp) WRITE(numout,*)'ice_rdgrft: ice ridging and rafting'
[8498]	137	IF(lwp) WRITE(numout,*)'~~~~~~~~~~'
	138	!
	139	zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE
	140	!
	141	ENDIF
[8426]	142
[8409]	143	!-----------------------------------------------------------------------------!
	144	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
	145	!-----------------------------------------------------------------------------!
	146	!
[8512]	147	CALL ice_rdgrft_prep ! prepare ridging
[8409]	148	!
[8486]	149	DO jj = 1, jpj ! Initialize arrays.
[8409]	150	DO ji = 1, jpi
	151	!-----------------------------------------------------------------------------!
	152	! 2) Dynamical inputs (closing rate, divu_adv, opning)
	153	!-----------------------------------------------------------------------------!
	154	!
	155	! 2.1 closing_net
	156	!-----------------
	157	! Compute the net rate of closing due to convergence
	158	! and shear, based on Flato and Hibler (1995).
	159	!
	160	! The energy dissipation rate is equal to the net closing rate
	161	! times the ice strength.
	162	!
	163	! NOTE: The NET closing rate is equal to the rate that open water
	164	! area is removed, plus the rate at which ice area is removed by
	165	! ridging, minus the rate at which area is added in new ridges.
	166	! The GROSS closing rate is equal to the first two terms (open
	167	! water closing and thin ice ridging) without the third term
	168	! (thick, newly ridged ice).
	169
[8515]	170	closing_net(ji,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
[8409]	171
	172	! 2.2 divu_adv
	173	!--------------
	174	! Compute divu_adv, the divergence rate given by the transport/
	175	! advection scheme, which may not be equal to divu as computed
	176	! from the velocity field.
	177	!
	178	! If divu_adv < 0, make sure the closing rate is large enough
	179	! to give asum = 1.0 after ridging.
	180
	181	divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
	182
	183	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
	184
	185	! 2.3 opning
	186	!------------
	187	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
	188	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
	189	END DO
	190	END DO
	191
	192	!-----------------------------------------------------------------------------!
	193	! 3) Ridging iteration
	194	!-----------------------------------------------------------------------------!
	195	niter = 1 ! iteration counter
	196	iterate_ridging = 1
	197
	198	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
	199
	200	! 3.2 closing_gross
	201	!-----------------------------------------------------------------------------!
	202	! Based on the ITD of ridging and ridged ice, convert the net
	203	! closing rate to a gross closing rate.
	204	! NOTE: 0 < aksum <= 1
	205	closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
	206
	207	! correction to closing rate and opening if closing rate is excessive
	208	!---------------------------------------------------------------------
	209	! Reduce the closing rate if more than 100% of the open water
	210	! would be removed. Reduce the opening rate proportionately.
	211	DO jj = 1, jpj
	212	DO ji = 1, jpi
	213	za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
	214	IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i
	215	opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice
	216	ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum
	217	opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice
	218	ENDIF
	219	END DO
	220	END DO
	221
	222	! correction to closing rate / opening if excessive ice removal
	223	!---------------------------------------------------------------
	224	! Reduce the closing rate if more than 100% of any ice category
	225	! would be removed. Reduce the opening rate proportionately.
	226	DO jl = 1, jpl
	227	DO jj = 1, jpj
	228	DO ji = 1, jpi
	229	za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
	230	IF( za > a_i(ji,jj,jl) ) THEN
[8498]	231	closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za
[8409]	232	ENDIF
	233	END DO
	234	END DO
	235	END DO
	236
	237	! 3.3 Redistribute area, volume, and energy.
	238	!-----------------------------------------------------------------------------!
	239	CALL ice_rdgrft_ridgeshift( opning, closing_gross )
	240
	241	! 3.4 Compute total area of ice plus open water after ridging.
	242	!-----------------------------------------------------------------------------!
	243	! This is in general not equal to one because of divergence during transport
	244	asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
	245
[8498]	246	! 3.5 Do we keep on iterating?
[8409]	247	!-----------------------------------------------------------------------------!
	248	! Check whether asum = 1. If not (because the closing and opening
	249	! rates were reduced above), ridge again with new rates.
	250	iterate_ridging = 0
	251	DO jj = 1, jpj
	252	DO ji = 1, jpi
	253	IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
	254	closing_net(ji,jj) = 0._wp
	255	opning (ji,jj) = 0._wp
	256	ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) )
	257	ELSE
	258	iterate_ridging = 1
	259	divu_adv (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
	260	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
	261	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
	262	ENDIF
	263	END DO
	264	END DO
	265	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
	266
	267	! Repeat if necessary.
	268	! NOTE: If strength smoothing is turned on, the ridging must be
[8486]	269	! iterated globally because of the boundary update in the smoothing.
[8409]	270	niter = niter + 1
[8486]	271	!
[8409]	272	IF( iterate_ridging == 1 ) THEN
[8512]	273	CALL ice_rdgrft_prep
[8409]	274	IF( niter > nitermax ) THEN
	275	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
	276	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
	277	ENDIF
	278	ENDIF
[8486]	279	!
[8409]	280	END DO !! on the do while over iter
	281
[8424]	282	CALL ice_var_agg( 1 )
[8409]	283
[8517]	284	! controls
	285	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
	286	IF( ln_ctl ) CALL ice_prt3D ('icerdgrft') ! prints
	287	IF( nn_timing == 1 ) CALL timing_stop ('icerdgrft') ! timing
[8409]	288	!
	289	END SUBROUTINE ice_rdgrft
	290
[8486]	291
[8512]	292	SUBROUTINE ice_rdgrft_prep
[8409]	293	!!---------------------------------------------------------------------!
[8512]	294	!! * ROUTINE ice_rdgrft_prep *
[8409]	295	!!
	296	!! ** Purpose : preparation for ridging and strength calculations
	297	!!
	298	!! ** Method : Compute the thickness distribution of the ice and open water
	299	!! participating in ridging and of the resulting ridges.
	300	!!---------------------------------------------------------------------!
[8486]	301	INTEGER :: ji, jj, jl ! dummy loop indices
[8498]	302	REAL(wp) :: z1_gstar, z1_astar, zhmean, zdummy ! local scalar
	303	REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n
[8409]	304	!------------------------------------------------------------------------------!
	305
[8498]	306	z1_gstar = 1._wp / rn_gstar
	307	z1_astar = 1._wp / rn_astar
[8409]	308
[8486]	309	CALL ice_var_zapsmall ! Zero out categories with very small areas
[8409]	310
[8498]	311	! ! Ice thickness needed for rafting
[8486]	312	WHERE( a_i(:,:,:) >= epsi20 ) ; ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)
	313	ELSEWHERE ; ht_i(:,:,:) = 0._wp
	314	END WHERE
[8409]	315
	316	!------------------------------------------------------------------------------!
	317	! 1) Participation function
	318	!------------------------------------------------------------------------------!
[8486]	319	!
[8409]	320	! Compute total area of ice plus open water.
	321	! This is in general not equal to one because of divergence during transport
	322	asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
[8486]	323	!
[8409]	324	! Compute cumulative thickness distribution function
[8498]	325	! Compute the cumulative thickness distribution function zGsum,
	326	! where zGsum(n) is the fractional area in categories 0 to n.
[8409]	327	! initial value (in h = 0) equals open water area
[8498]	328	zGsum(:,:,-1) = 0._wp
	329	zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:)
[8409]	330	DO jl = 1, jpl
[8498]	331	zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl)
[8409]	332	END DO
	333
	334	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
	335	!--------------------------------------------------------------------------------------------------
	336	! Compute the participation function athorn; this is analogous to
	337	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
	338	! area lost from category n due to ridging/closing
	339	! athorn(n) = total area lost due to ridging/closing
	340	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
	341	!
	342	! The expressions for athorn are found by integrating b(h)g(h) between
	343	! the category boundaries.
	344	! athorn is always >= 0 and SUM(athorn(0:jpl))=1
	345	!-----------------------------------------------------------------
[8486]	346	!
[8512]	347	IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975)
[8409]	348	DO jl = 0, jpl
	349	DO jj = 1, jpj
	350	DO ji = 1, jpi
[8498]	351	IF ( zGsum(ji,jj,jl) < rn_gstar ) THEN
	352	athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
	353	& ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
	354	ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
	355	athorn(ji,jj,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jj,jl-1) ) * &
	356	& ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar ) * z1_gstar )
[8409]	357	ELSE
	358	athorn(ji,jj,jl) = 0._wp
	359	ENDIF
	360	END DO
	361	END DO
	362	END DO
[8486]	363	!
[8512]	364	ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007)
[8409]	365	!
[8498]	366	zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) ) ! precompute exponential terms using zGsum as a work array
[8409]	367	DO jl = -1, jpl
[8498]	368	zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy
[8409]	369	END DO
	370	DO jl = 0, jpl
[8498]	371	athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
[8409]	372	END DO
	373	!
	374	ENDIF
	375
[8486]	376	! !--- Ridging and rafting participation concentrations
	377	IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting
[8409]	378	DO jl = 1, jpl
	379	DO jj = 1, jpj
	380	DO ji = 1, jpi
[8498]	381	aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl)
[8409]	382	araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
	383	END DO
	384	END DO
	385	END DO
[8498]	386	ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone
	387	aridge(:,:,:) = athorn(:,:,1:jpl)
	388	araft (:,:,:) = 0._wp
	389	ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone
	390	aridge(:,:,:) = 0._wp
	391	araft (:,:,:) = athorn(:,:,1:jpl)
	392	ELSE !- no ridging & no rafting
	393	aridge(:,:,:) = 0._wp
	394	araft (:,:,:) = 0._wp
[8409]	395	ENDIF
	396
	397	!-----------------------------------------------------------------
	398	! 2) Transfer function
	399	!-----------------------------------------------------------------
	400	! Compute max and min ridged ice thickness for each ridging category.
	401	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
	402	!
	403	! This parameterization is a modified version of Hibler (1980).
[8498]	404	! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
[8409]	405	! and for very thick ridging ice must be >= krdgmin*hi
	406	!
	407	! The minimum ridging thickness, hrmin, is equal to 2*hi
	408	! (i.e., rafting) and for very thick ridging ice is
[8498]	409	! constrained by hrmin <= (zhmean + hi)/2.
[8409]	410	!
	411	! The maximum ridging thickness, hrmax, is determined by
[8498]	412	! zhmean and hrmin.
[8409]	413	!
	414	! These modifications have the effect of reducing the ice strength
	415	! (relative to the Hibler formulation) when very thick ice is
	416	! ridging.
	417	!
	418	! aksum = net area removed/ total area removed
	419	! where total area removed = area of ice that ridges
	420	! net area removed = total area removed - area of new ridges
	421	!-----------------------------------------------------------------
	422
	423	aksum(:,:) = athorn(:,:,0)
	424	! Transfer function
	425	DO jl = 1, jpl !all categories have a specific transfer function
	426	DO jj = 1, jpj
	427	DO ji = 1, jpi
[8486]	428	IF ( athorn(ji,jj,jl) > 0._wp ) THEN
[8498]	429	zhmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
	430	hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
	431	hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
[8409]	432	hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
[8498]	433	krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
[8486]	434	!
[8409]	435	! Normalization factor : aksum, ensures mass conservation
	436	aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) ) &
	437	& + araft (ji,jj,jl) * ( 1._wp - kraft )
	438	ELSE
	439	hrmin(ji,jj,jl) = 0._wp
	440	hrmax(ji,jj,jl) = 0._wp
	441	hraft(ji,jj,jl) = 0._wp
	442	krdg (ji,jj,jl) = 1._wp
	443	ENDIF
	444	END DO
	445	END DO
	446	END DO
	447	!
[8512]	448	END SUBROUTINE ice_rdgrft_prep
[8409]	449
	450
	451	SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross )
	452	!!----------------------------------------------------------------------
[8512]	453	!! * ROUTINE ice_rdgrft_strength *
[8409]	454	!!
	455	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
	456	!!
	457	!! ** Method : Remove area, volume, and energy from each ridging category
	458	!! and add to thicker ice categories.
	459	!!----------------------------------------------------------------------
	460	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
[8498]	461	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area retreats, excluding area of new ridges
[8409]	462	!
	463	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
	464	INTEGER :: ij ! horizontal index, combines i and j loops
	465	INTEGER :: icells ! number of cells with a_i > puny
[8498]	466	REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2
[8409]	467
	468	INTEGER , DIMENSION(jpij) :: indxi, indxj ! compressed indices
[8498]	469	REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2
[8409]	470
	471	REAL(wp), DIMENSION(jpij) :: afrac ! fraction of category area ridged
	472	REAL(wp), DIMENSION(jpij) :: ardg1 , ardg2 ! area of ice ridged & new ridges
	473	REAL(wp), DIMENSION(jpij) :: vsrdg , esrdg ! snow volume & energy of ridging ice
	474	! MV MP 2016
	475	REAL(wp), DIMENSION(jpij) :: vprdg ! pond volume of ridging ice
	476	REAL(wp), DIMENSION(jpij) :: aprdg1 ! pond area of ridging ice
	477	REAL(wp), DIMENSION(jpij) :: aprdg2 ! pond area of ridging ice
	478	! END MV MP 2016
	479	REAL(wp), DIMENSION(jpij) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2
	480
	481	REAL(wp), DIMENSION(jpij) :: vrdg1 ! volume of ice ridged
	482	REAL(wp), DIMENSION(jpij) :: vrdg2 ! volume of new ridges
	483	REAL(wp), DIMENSION(jpij) :: vsw ! volume of seawater trapped into ridges
	484	REAL(wp), DIMENSION(jpij) :: srdg1 ! sal*volume of ice ridged
	485	REAL(wp), DIMENSION(jpij) :: srdg2 ! sal*volume of new ridges
	486	REAL(wp), DIMENSION(jpij) :: smsw ! sal*volume of water trapped into ridges
	487	REAL(wp), DIMENSION(jpij) :: oirdg1, oirdg2 ! ice age of ice ridged
	488
	489	REAL(wp), DIMENSION(jpij) :: afrft ! fraction of category area rafted
	490	REAL(wp), DIMENSION(jpij) :: arft1 , arft2 ! area of ice rafted and new rafted zone
	491	REAL(wp), DIMENSION(jpij) :: virft , vsrft ! ice & snow volume of rafting ice
	492	! MV MP 2016
	493	REAL(wp), DIMENSION(jpij) :: vprft ! pond volume of rafting ice
	494	REAL(wp), DIMENSION(jpij) :: aprft1 ! pond area of rafted ice
	495	REAL(wp), DIMENSION(jpij) :: aprft2 ! pond area of new rafted ice
	496	! END MV MP 2016
	497	REAL(wp), DIMENSION(jpij) :: esrft , smrft ! snow energy & salinity of rafting ice
	498	REAL(wp), DIMENSION(jpij) :: oirft1, oirft2 ! ice age of ice rafted
	499
	500	REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice
	501	REAL(wp), DIMENSION(jpij,nlay_i) :: erdg1 ! enth*volume of ice ridged
	502	REAL(wp), DIMENSION(jpij,nlay_i) :: erdg2 ! enth*volume of new ridges
	503	REAL(wp), DIMENSION(jpij,nlay_i) :: ersw ! enth of water trapped into ridges
	504	!!----------------------------------------------------------------------
	505
	506	!-------------------------------------------------------------------------------
	507	! 1) Compute change in open water area due to closing and opening.
	508	!-------------------------------------------------------------------------------
	509	DO jj = 1, jpj
	510	DO ji = 1, jpi
	511	ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) + &
	512	& ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
	513	END DO
	514	END DO
	515
	516	!-----------------------------------------------------------------
	517	! 3) Pump everything from ice which is being ridged / rafted
	518	!-----------------------------------------------------------------
	519	! Compute the area, volume, and energy of ice ridging in each
	520	! category, along with the area of the resulting ridge.
	521
	522	DO jl1 = 1, jpl !jl1 describes the ridging category
	523
	524	!------------------------------------------------
	525	! 3.1) Identify grid cells with nonzero ridging
	526	!------------------------------------------------
	527	icells = 0
	528	DO jj = 1, jpj
	529	DO ji = 1, jpi
	530	IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
	531	icells = icells + 1
	532	indxi(icells) = ji
	533	indxj(icells) = jj
	534	ENDIF
	535	END DO
	536	END DO
	537
	538	DO ij = 1, icells
	539	ji = indxi(ij) ; jj = indxj(ij)
	540
	541	!--------------------------------------------------------------------
	542	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
	543	!--------------------------------------------------------------------
	544	ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
	545	arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
	546
	547	!---------------------------------------------------------------
	548	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
	549	!---------------------------------------------------------------
	550	afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
	551	afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
	552	ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
	553	arft2(ij) = arft1(ij) * kraft
	554
	555	!--------------------------------------------------------------------------
	556	! 3.4) Substract area, volume, and energy from ridging
	557	! / rafting category n1.
	558	!--------------------------------------------------------------------------
	559	vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
[8512]	560	vrdg2(ij) = vrdg1(ij) * ( 1. + rn_porordg )
	561	vsw (ij) = vrdg1(ij) * rn_porordg
[8409]	562
	563	vsrdg (ij) = v_s (ji,jj, jl1) * afrac(ij)
	564	esrdg (ij) = e_s (ji,jj,1,jl1) * afrac(ij)
	565	!MV MP 2016
	566	IF ( nn_pnd_scheme > 0 ) THEN
	567	aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij)
	568	aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
	569	vprdg(ij) = v_ip(ji,jj, jl1) * afrac(ij)
	570	ENDIF
	571	! END MV MP 2016
	572	srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij)
	573	oirdg1(ij) = oa_i (ji,jj, jl1) * afrac(ij)
	574	oirdg2(ij) = oa_i (ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
	575
	576	! rafting volumes, heat contents ...
	577	virft (ij) = v_i (ji,jj, jl1) * afrft(ij)
	578	vsrft (ij) = v_s (ji,jj, jl1) * afrft(ij)
	579	!MV MP 2016
	580	IF ( nn_pnd_scheme > 0 ) THEN
	581	aprft1(ij) = a_ip (ji,jj, jl1) * afrft(ij)
	582	aprft2(ij) = a_ip (ji,jj, jl1) * afrft(ij) * kraft
	583	vprft(ij) = v_ip(ji,jj,jl1) * afrft(ij)
	584	ENDIF
	585	! END MV MP 2016
	586	srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij)
	587	esrft (ij) = e_s (ji,jj,1,jl1) * afrft(ij)
	588	smrft (ij) = smv_i(ji,jj, jl1) * afrft(ij)
	589	oirft1(ij) = oa_i (ji,jj, jl1) * afrft(ij)
	590	oirft2(ij) = oa_i (ji,jj, jl1) * afrft(ij) * kraft
	591
	592	!-----------------------------------------------------------------
	593	! 3.5) Compute properties of new ridges
	594	!-----------------------------------------------------------------
	595	smsw(ij) = vsw(ij) * sss_m(ji,jj) ! salt content of seawater frozen in voids
	596	srdg2(ij) = srdg1(ij) + smsw(ij) ! salt content of new ridge
	597
	598	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
	599	wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids
	600
	601	! virtual salt flux to keep salinity constant
	602	IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN
	603	srdg2(ij) = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) ) ! ridge salinity = sm_i
	604	sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj) * vsw(ij) * rhoic * r1_rdtice & ! put back sss_m into the ocean
	605	& - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean
	606	ENDIF
	607
	608	!------------------------------------------
	609	! 3.7 Put the snow somewhere in the ocean
	610	!------------------------------------------
	611	! Place part of the snow lost by ridging into the ocean.
	612	! Note that esrdg > 0; the ocean must cool to melt snow.
	613	! If the ocean temp = Tf already, new ice must grow.
	614	! During the next time step, thermo_rates will determine whether
	615	! the ocean cools or new ice grows.
[8512]	616	wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnwrdg ) &
	617	& + rhosn * vsrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! fresh water source for ocean
[8409]	618
[8512]	619	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnwrdg ) &
	620	& - esrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
[8409]	621
	622	! MV MP 2016
	623	!------------------------------------------
	624	! 3.X Put the melt pond water in the ocean
	625	!------------------------------------------
	626	! Place part of the melt pond volume into the ocean.
	627	IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
[8512]	628	wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpndrdg ) &
	629	& + rhofw * vprft(ij) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice ! fresh water source for ocean
[8409]	630	ENDIF
	631	! END MV MP 2016
	632
	633	!-----------------------------------------------------------------
	634	! 3.8 Compute quantities used to apportion ice among categories
	635	! in the n2 loop below
	636	!-----------------------------------------------------------------
	637	dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) )
	638	dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
	639
	640
	641	! update jl1 (removing ridged/rafted area)
	642	a_i (ji,jj, jl1) = a_i (ji,jj, jl1) - ardg1 (ij) - arft1 (ij)
	643	v_i (ji,jj, jl1) = v_i (ji,jj, jl1) - vrdg1 (ij) - virft (ij)
	644	v_s (ji,jj, jl1) = v_s (ji,jj, jl1) - vsrdg (ij) - vsrft (ij)
	645	e_s (ji,jj,1,jl1) = e_s (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
	646	smv_i(ji,jj, jl1) = smv_i(ji,jj, jl1) - srdg1 (ij) - smrft (ij)
	647	oa_i (ji,jj, jl1) = oa_i (ji,jj, jl1) - oirdg1(ij) - oirft1(ij)
	648
	649	! MV MP 2016
	650	IF ( nn_pnd_scheme > 0 ) THEN
	651	v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
	652	a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij)
	653	ENDIF
	654	! END MV MP 2016
	655
	656	END DO
	657
	658	!--------------------------------------------------------------------
	659	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
	660	! compute ridged ice enthalpy
	661	!--------------------------------------------------------------------
	662	DO jk = 1, nlay_i
	663	DO ij = 1, icells
	664	ji = indxi(ij) ; jj = indxj(ij)
	665	! heat content of ridged ice
	666	erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij)
	667	eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)
	668
	669	! enthalpy of the trapped seawater (J/m2, >0)
	670	! clem: if sst>0, then ersw <0 (is that possible?)
	671	ersw(ij,jk) = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
	672
	673	! heat flux to the ocean
	674	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
	675
	676	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
	677	erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
	678
	679	! update jl1
[8498]	680	e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
[8409]	681
	682	END DO
	683	END DO
	684
	685	!-------------------------------------------------------------------------------
	686	! 4) Add area, volume, and energy of new ridge to each category jl2
	687	!-------------------------------------------------------------------------------
	688	DO jl2 = 1, jpl
	689	! over categories to which ridged/rafted ice is transferred
	690	DO ij = 1, icells
	691	ji = indxi(ij) ; jj = indxj(ij)
	692
	693	! Compute the fraction of ridged ice area and volume going to thickness category jl2.
	694	IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
	695	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
	696	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
	697	farea = ( hR - hL ) * dhr(ij)
	698	fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
	699	ELSE
	700	farea = 0._wp
	701	fvol(ij) = 0._wp
	702	ENDIF
	703
	704	! Compute the fraction of rafted ice area and volume going to thickness category jl2
[8486]	705	!!gm see above IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
	706	IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2) ) THEN ; zswitch(ij) = 1._wp
	707	ELSE ; zswitch(ij) = 0._wp
[8409]	708	ENDIF
[8486]	709	!
[8409]	710	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ( ardg2 (ij) * farea + arft2 (ij) * zswitch(ij) )
	711	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + ( oirdg2(ij) * farea + oirft2(ij) * zswitch(ij) )
	712	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
	713	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
[8512]	714	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + ( vsrdg (ij) * rn_fsnwrdg * fvol(ij) + &
	715	& vsrft (ij) * rn_fsnwrft * zswitch(ij) )
	716	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnwrdg * fvol(ij) + &
	717	& esrft (ij) * rn_fsnwrft * zswitch(ij) )
[8409]	718	! MV MP 2016
	719	IF ( nn_pnd_scheme > 0 ) THEN
[8512]	720	v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + ( vprdg (ij) * rn_fpndrdg * fvol (ij) &
	721	& + vprft (ij) * rn_fpndrft * zswitch(ij) )
	722	a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + ( aprdg2(ij) * rn_fpndrdg * farea &
	723	& + aprft2(ij) * rn_fpndrft * zswitch(ji) )
[8409]	724	ENDIF
	725	! END MV MP 2016
	726	END DO
	727
	728	! Transfer ice energy to category jl2 by ridging
	729	DO jk = 1, nlay_i
	730	DO ij = 1, icells
	731	ji = indxi(ij) ; jj = indxj(ij)
	732	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)
	733	END DO
	734	END DO
	735	!
	736	END DO ! jl2
[8486]	737	!
[8409]	738	END DO ! jl1 (deforming categories)
	739	!
	740	END SUBROUTINE ice_rdgrft_ridgeshift
	741
[8486]	742
[8512]	743	SUBROUTINE ice_rdgrft_strength
[8409]	744	!!----------------------------------------------------------------------
[8512]	745	!! * ROUTINE ice_rdgrft_strength *
[8409]	746	!!
	747	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
	748	!!
	749	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
	750	!! dissipated per unit area removed from the ice pack under compression,
	751	!! and assumed proportional to the change in potential energy caused
	752	!! by ridging. Note that only Hibler's formulation is stable and that
	753	!! ice strength has to be smoothed
	754	!!----------------------------------------------------------------------
	755	INTEGER :: ji,jj, jl ! dummy loop indices
[8498]	756	INTEGER :: ismooth ! smoothing the resistance to deformation
	757	INTEGER :: itframe ! number of time steps for the P smoothing
[8409]	758	REAL(wp) :: zp, z1_3 ! local scalars
	759	REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here
	760	REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps
	761	!!----------------------------------------------------------------------
	762
[8498]	763	! !--------------------------------------------------!
[8512]	764	CALL ice_rdgrft_prep ! Thickness distribution of ridging and ridged ice !
[8498]	765	! !--------------------------------------------------!
[8409]	766
[8498]	767	! !--------------------------------------------------!
[8514]	768	IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method !
[8498]	769	! !--------------------------------------------------!
[8409]	770	z1_3 = 1._wp / 3._wp
	771	DO jl = 1, jpl
[8498]	772	WHERE( athorn(:,:,jl) > 0._wp )
	773	strength(:,:) = - athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl) & ! PE loss from deforming ice
	774	& + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl) & ! PE gain from rafting ice
	775	& + aridge(:,:,jl) * krdg(:,:,jl) * z1_3 * & ! PE gain from ridging ice
	776	& ( hrmax(:,:,jl) * hrmax(:,:,jl) + &
	777	& hrmin(:,:,jl) * hrmin(:,:,jl) + &
	778	& hrmax(:,:,jl) * hrmin(:,:,jl) )
	779	ELSEWHERE
	780	strength(:,:) = 0._wp
	781	END WHERE
[8409]	782	END DO
[8512]	783	strength(:,:) = rn_perdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)
	784	! where zdrho = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_perdg accounts for frictional dissipation
[8498]	785	ismooth = 1
	786	! !--------------------------------------------------!
[8514]	787	ELSEIF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method !
[8498]	788	! !--------------------------------------------------!
[8517]	789	strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
[8409]	790	!
[8498]	791	ismooth = 1
[8409]	792	!
[8498]	793	ENDIF
	794	! !--------------------------------------------------!
	795	SELECT CASE( ismooth ) ! Smoothing ice strength !
	796	! !--------------------------------------------------!
	797	CASE( 1 ) !--- Spatial smoothing
[8409]	798	DO jj = 2, jpjm1
	799	DO ji = 2, jpim1
	800	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
	801	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
	802	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
	803	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
	804	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
	805	ELSE
	806	zworka(ji,jj) = 0._wp
	807	ENDIF
	808	END DO
	809	END DO
[8498]	810
[8409]	811	DO jj = 2, jpjm1
	812	DO ji = 2, jpim1
	813	strength(ji,jj) = zworka(ji,jj)
	814	END DO
	815	END DO
	816	CALL lbc_lnk( strength, 'T', 1. )
[8486]	817	!
[8498]	818	CASE( 2 ) !--- Temporal smoothing
[8409]	819	IF ( kt_ice == nit000 ) THEN
	820	zstrp1(:,:) = 0._wp
	821	zstrp2(:,:) = 0._wp
	822	ENDIF
[8486]	823	!
[8409]	824	DO jj = 2, jpjm1
	825	DO ji = 2, jpim1
	826	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN
[8498]	827	itframe = 1 ! number of time steps for the running mean
	828	IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
	829	IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
	830	zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
[8486]	831	zstrp2 (ji,jj) = zstrp1 (ji,jj)
	832	zstrp1 (ji,jj) = strength(ji,jj)
[8409]	833	strength(ji,jj) = zp
	834	ENDIF
	835	END DO
	836	END DO
[8498]	837	CALL lbc_lnk( strength, 'T', 1. )
[8486]	838	!
	839	END SELECT
[8409]	840	!
[8512]	841	END SUBROUTINE ice_rdgrft_strength
[8409]	842
[8486]	843
[8409]	844	SUBROUTINE ice_rdgrft_init
	845	!!-------------------------------------------------------------------
[8486]	846	!! * ROUTINE ice_rdgrft_init *
[8409]	847	!!
	848	!! ** Purpose : Physical constants and parameters linked
	849	!! to the mechanical ice redistribution
	850	!!
[8512]	851	!! ** Method : Read the namice_rdgrft namelist
[8409]	852	!! and check the parameters values
	853	!! called at the first timestep (nit000)
	854	!!
[8512]	855	!! ** input : Namelist namice_rdgrft
[8409]	856	!!-------------------------------------------------------------------
	857	INTEGER :: ios ! Local integer output status for namelist read
[8486]	858	!!
[8514]	859	NAMELIST/namice_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
	860	& ln_str_R75, rn_perdg, &
[8515]	861	& rn_csrdg , &
[8512]	862	& ln_partf_lin, rn_gstar, &
	863	& ln_partf_exp, rn_astar, &
	864	& ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, &
	865	& ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft
[8409]	866	!!-------------------------------------------------------------------
	867	!
	868	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
[8512]	869	READ ( numnam_ice_ref, namice_rdgrft, IOSTAT = ios, ERR = 901)
	870	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in reference namelist', lwp )
[8486]	871	!
[8512]	872	REWIND( numnam_ice_cfg ) ! Namelist namice_rdgrft in configuration namelist : Ice mechanical ice redistribution
	873	READ ( numnam_ice_cfg, namice_rdgrft, IOSTAT = ios, ERR = 902 )
	874	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in configuration namelist', lwp )
	875	IF(lwm) WRITE ( numoni, namice_rdgrft )
[8409]	876	!
	877	IF (lwp) THEN ! control print
	878	WRITE(numout,*)
[8514]	879	WRITE(numout,*) 'ice_rdgrft_init: ice parameters for ridging/rafting '
[8498]	880	WRITE(numout,*) '~~~~~~~~~~~~~~~'
[8512]	881	WRITE(numout,*) ' Namelist namice_rdgrft'
[8514]	882	WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79
[8512]	883	WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar
	884	WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg
[8514]	885	WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75
[8512]	886	WRITE(numout,*) ' Ratio of ridging work to PotEner change in ridging rn_perdg = ', rn_perdg
[8515]	887	WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg
[8512]	888	WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin
	889	WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar
	890	WRITE(numout,*) ' Exponential ridging participation function ln_partf_exp = ', ln_partf_exp
	891	WRITE(numout,) ' Equivalent to G for an exponential function rn_astar = ', rn_astar
	892	WRITE(numout,*) ' Ridging of ice sheets or not ln_ridging = ', ln_ridging
	893	WRITE(numout,*) ' max ridged ice thickness rn_hstar = ', rn_hstar
	894	WRITE(numout,*) ' Initial porosity of ridges rn_porordg = ', rn_porordg
	895	WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg
	896	WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg
	897	WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting
	898	WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
	899	WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
	900	WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft
	901	WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft
[8409]	902	ENDIF
	903	!
[8514]	904	IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN
	905	CALL ctl_stop( 'ice_rdgrft_init: choose one and only one formulation for ice strength (ln_str_H79 or ln_str_R75)' )
[8512]	906	ENDIF
	907	!
	908	IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
	909	CALL ctl_stop( 'ice_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
	910	ENDIF
	911	!
	912	END SUBROUTINE ice_rdgrft_init
[8409]	913
	914	#else
	915	!!----------------------------------------------------------------------
	916	!! Default option Empty module NO LIM-3 sea-ice model
	917	!!----------------------------------------------------------------------
	918	#endif
[8486]	919
[8409]	920	!!======================================================================
	921	END MODULE icerdgrft

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

Download in other formats: