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.

iceitd.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/iceitd.F90 @ 8486

Last change on this file since 8486 was 8486, checked in by clem, 7 years ago
changes in style - part1 - (now the code looks better txs to Gurvan's comments)
File size: 31.1 KB

Rev	Line
[8422]	1	MODULE iceitd
	2	!!======================================================================
	3	!! * MODULE iceitd *
	4	!! LIM3 ice model : ice thickness distribution: Thermodynamics
	5	!!======================================================================
	6	!! History : - ! (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
	7	!! 3.0 ! 2005-12 (M. Vancoppenolle) adaptation to LIM-3
	8	!! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age
	9	!! - ! 2007-04 (M. Vancoppenolle) Mass conservation checked
	10	!!----------------------------------------------------------------------
	11	#if defined key_lim3
	12	!!----------------------------------------------------------------------
[8486]	13	!! 'key_lim3' LIM3 sea-ice model
[8422]	14	!!----------------------------------------------------------------------
	15	!! ice_itd_rem :
	16	!! ice_itd_reb :
	17	!! ice_itd_glinear :
	18	!! ice_itd_shiftice :
	19	!!----------------------------------------------------------------------
[8486]	20	USE par_oce ! ocean parameters
	21	USE dom_oce ! ocean domain
	22	USE phycst ! physical constants
	23	USE ice1D ! sea-ice: thermodynamic variables
	24	USE ice ! sea-ice: variables
	25	USE icectl ! sea-ice: conservation tests
	26	USE icetab ! sea-ice: convert 1D<=>2D
[8422]	27	!
[8486]	28	USE prtctl ! Print control
	29	USE in_out_manager ! I/O manager
	30	USE lib_mpp ! MPP library
	31	USE lib_fortran ! to use key_nosignedzero
[8422]	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC ice_itd_rem ! called in icethd
	37	PUBLIC ice_itd_reb ! called in iceerr
	38
	39	!!----------------------------------------------------------------------
[8486]	40	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8422]	41	!! $Id: iceitd.F90 8420 2017-08-08 12:18:46Z clem $
	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
	46	SUBROUTINE ice_itd_rem( kt )
	47	!!------------------------------------------------------------------
	48	!! * ROUTINE ice_itd_rem *
	49	!!
	50	!! ** Purpose : computes the redistribution of ice thickness
	51	!! after thermodynamic growth of ice thickness
	52	!!
	53	!! ** Method : Linear remapping
	54	!!
	55	!! References : W.H. Lipscomb, JGR 2001
	56	!!------------------------------------------------------------------
	57	INTEGER , INTENT (in) :: kt ! Ocean time step
	58	!
	59	INTEGER :: ji, jj, jl, jcat ! dummy loop index
	60	INTEGER :: nidx2 ! local integer
	61	REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars
	62	REAL(wp) :: zx2, zwk2, zda0, zetamax ! - -
	63	REAL(wp) :: zx3
	64	REAL(wp) :: zslope ! used to compute local thermodynamic "speeds"
	65	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
	66
	67	INTEGER , DIMENSION(jpij) :: idxice2 ! compute remapping or not
	68	INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index
	69	REAL(wp), DIMENSION(jpij,jpl) :: zdhice ! ice thickness increment
	70	REAL(wp), DIMENSION(jpij,jpl) :: g0, g1 ! coefficients for fitting the line of the ITD
	71	REAL(wp), DIMENSION(jpij,jpl) :: hL, hR ! left and right boundary for the ITD for each thickness
	72	REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! local increment of ice area and volume
	73	REAL(wp), DIMENSION(jpij) :: zhb0, zhb1 ! category boundaries for thinnes categories
	74	REAL(wp), DIMENSION(jpij,0:jpl) :: zhbnew ! new boundaries of ice categories
	75	!!------------------------------------------------------------------
	76
	77	IF( kt == nit000 .AND. lwp) THEN
	78	WRITE(numout,*)
	79	WRITE(numout,*) 'ice_itd_rem : Remapping the ice thickness distribution'
	80	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	81	ENDIF
	82
	83	IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	84
	85	!-----------------------------------------------------------------------------------------------
	86	! 1) Identify grid cells with ice
	87	!-----------------------------------------------------------------------------------------------
[8486]	88	nidx = 0 ; idxice(:) = 0
[8422]	89	DO jj = 1, jpj
	90	DO ji = 1, jpi
	91	IF ( at_i(ji,jj) > epsi10 ) THEN
	92	nidx = nidx + 1
	93	idxice( nidx ) = (jj - 1) * jpi + ji
	94	ENDIF
	95	END DO
	96	END DO
	97
	98	!-----------------------------------------------------------------------------------------------
	99	! 2) Compute new category boundaries
	100	!-----------------------------------------------------------------------------------------------
	101	IF( nidx > 0 ) THEN
	102
	103	zdhice(:,:) = 0._wp
	104	zhbnew(:,:) = 0._wp
	105
	106	CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i )
	107	CALL tab_3d_2d( nidx, idxice(1:nidx), ht_ib_2d(1:nidx,1:jpl), ht_i_b )
	108	CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )
	109	CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b )
	110
	111	DO jl = 1, jpl
	112	! Compute thickness change in each ice category
	113	DO ji = 1, nidx
	114	zdhice(ji,jl) = ht_i_2d(ji,jl) - ht_ib_2d(ji,jl)
	115	END DO
	116	END DO
	117
	118	! --- New boundaries for category 1:jpl-1 --- !
	119	DO jl = 1, jpl - 1
[8486]	120	!
[8422]	121	DO ji = 1, nidx
	122	!
	123	! --- New boundary: Hn* = Hn + Fn*dt --- !
	124	! Fndt = ( fn + (fn+1 - fn)/(hn+1 - hn) (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b)
	125	!
	126	IF ( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl+1) & a(jl) /= 0
	127	zslope = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( ht_ib_2d(ji,jl+1) - ht_ib_2d(ji,jl) )
	128	zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - ht_ib_2d(ji,jl) )
	129	ELSEIF( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN ! a(jl+1)=0 => Hn* = Hn + fn*dt
	130	zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
	131	ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl)=0 => Hn* = Hn + fn+1*dt
	132	zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
	133	ELSE ! a(jl+1) & a(jl) = 0
	134	zhbnew(ji,jl) = hi_max(jl)
	135	ENDIF
[8486]	136	!
[8422]	137	! --- 2 conditions for remapping --- !
	138	! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi
	139	! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
	140	! in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
	141	IF( a_i_2d(ji,jl ) > epsi10 .AND. ht_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) idxice(ji) = 0
	142	IF( a_i_2d(ji,jl+1) > epsi10 .AND. ht_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) idxice(ji) = 0
	143
	144	! 2) Hn-1 < Hn* < Hn+1
	145	IF( zhbnew(ji,jl) < hi_max(jl-1) ) idxice(ji) = 0
	146	IF( zhbnew(ji,jl) > hi_max(jl+1) ) idxice(ji) = 0
	147
	148	END DO
	149	END DO
[8486]	150	!
[8422]	151	! --- New boundaries for category jpl --- !
	152	DO ji = 1, nidx
	153	IF( a_i_2d(ji,jpl) > epsi10 ) THEN
	154	zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
	155	ELSE
	156	zhbnew(ji,jpl) = hi_max(jpl)
	157	ENDIF
	158
	159	! --- 1 additional condition for remapping (1st category) --- !
	160	! H0+epsi < h1(t) < H1-epsi
	161	! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
	162	! in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
	163	IF( ht_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) idxice(ji) = 0
	164	IF( ht_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) idxice(ji) = 0
	165	END DO
[8486]	166	!
[8422]	167	!-----------------------------------------------------------------------------------------------
	168	! 3) Identify cells where remapping
	169	!-----------------------------------------------------------------------------------------------
	170	nidx2 = 0 ; idxice2(:) = 0
	171	DO ji = 1, nidx
	172	IF( idxice(ji) /= 0 ) THEN
	173	nidx2 = nidx2 + 1
	174	idxice2(nidx2) = idxice(ji)
	175	zhbnew(nidx2,:) = zhbnew(ji,:) ! adjust zhbnew to new indices
	176	ENDIF
	177	END DO
	178	idxice(:) = idxice2(:)
	179	nidx = nidx2
[8486]	180	!
[8422]	181	ENDIF
	182
	183
	184	!-----------------------------------------------------------------------------------------------
	185	! 4) Compute g(h)
	186	!-----------------------------------------------------------------------------------------------
	187	IF( nidx > 0 ) THEN
[8486]	188	!
	189	zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1)
	190	g0(:,:) = 0._wp ; g1(:,:) = 0._wp
	191	hL(:,:) = 0._wp ; hR(:,:) = 0._wp
	192	!
[8422]	193	DO jl = 1, jpl
[8486]	194	!
[8422]	195	CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) )
	196	CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) )
	197	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )
	198	CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )
[8486]	199	!
[8422]	200	IF( jl == 1 ) THEN
[8486]	201	!
[8422]	202	! --- g(h) for category 1 --- !
[8486]	203	CALL ice_itd_glinear( zhb0(1:nidx) , zhb1(1:nidx) , ht_ib_1d(1:nidx) , a_i_1d(1:nidx) , & ! in
	204	& g0 (1:nidx,1), g1 (1:nidx,1), hL (1:nidx,1), hR (1:nidx,1) ) ! out
	205	!
[8422]	206	! Area lost due to melting of thin ice
	207	DO ji = 1, nidx
[8486]	208	!
[8422]	209	IF( a_i_1d(ji) > epsi10 ) THEN
[8486]	210	!
[8422]	211	zdh0 = ht_i_1d(ji) - ht_ib_1d(ji)
	212	IF( zdh0 < 0.0 ) THEN !remove area from category 1
	213	zdh0 = MIN( -zdh0, hi_max(1) )
	214	!Integrate g(1) from 0 to dh0 to estimate area melted
	215	zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
[8486]	216	!
[8422]	217	IF( zetamax > 0.0 ) THEN
	218	zx1 = zetamax
	219	zx2 = 0.5 * zetamax * zetamax
	220	zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed
	221	zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i
	222	zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting
	223	! of thin ice (zdamax > 0)
	224	! Remove area, conserving volume
	225	ht_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
	226	a_i_1d(ji) = a_i_1d(ji) - zda0
	227	v_i_1d(ji) = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ?
	228	ENDIF
[8486]	229	!
[8422]	230	ELSE ! if ice accretion zdh0 > 0
	231	! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
	232	zhbnew(ji,0) = MIN( zdh0, hi_max(1) )
	233	ENDIF
[8486]	234	!
[8422]	235	ENDIF
[8486]	236	!
[8422]	237	END DO
[8486]	238	!
[8422]	239	CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) )
	240	CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )
	241	CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )
[8486]	242	!
[8422]	243	ENDIF ! jl=1
[8486]	244	!
[8422]	245	! --- g(h) for each thickness category --- !
[8486]	246	CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx) , a_i_1d(1:nidx) , & ! in
	247	& g0 (1:nidx,jl ), g1 (1:nidx,jl), hL (1:nidx,jl), hR (1:nidx,jl) ) ! out
	248	!
[8422]	249	END DO
	250
	251	!-----------------------------------------------------------------------------------------------
	252	! 5) Compute area and volume to be shifted across each boundary (Eq. 18)
	253	!-----------------------------------------------------------------------------------------------
	254	DO jl = 1, jpl - 1
[8486]	255	!
[8422]	256	DO ji = 1, nidx
[8486]	257	!
[8422]	258	! left and right integration limits in eta space
	259	IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
	260	zetamin = MAX( hi_max(jl) , hL(ji,jl) ) - hL(ji,jl) ! hi_max(jl) - hL
	261	zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl) ! hR - hL
	262	jdonor(ji,jl) = jl
	263	ELSE ! Hn* <= Hn => transfer from jl+1 to jl
	264	zetamin = 0.0
	265	zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1) ! hi_max(jl) - hL
	266	jdonor(ji,jl) = jl + 1
	267	ENDIF
	268	zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
[8486]	269	!
[8422]	270	zx1 = zetamax - zetamin
	271	zwk1 = zetamin * zetamin
	272	zwk2 = zetamax * zetamax
	273	zx2 = 0.5 * ( zwk2 - zwk1 )
	274	zwk1 = zwk1 * zetamin
	275	zwk2 = zwk2 * zetamax
	276	zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 )
	277	jcat = jdonor(ji,jl)
	278	zdaice(ji,jl) = g1(ji,jcat)zx2 + g0(ji,jcat)zx1
	279	zdvice(ji,jl) = g1(ji,jcat)zx3 + g0(ji,jcat)zx2 + zdaice(ji,jl)*hL(ji,jcat)
[8486]	280	!
[8422]	281	END DO
	282	END DO
	283
	284	!----------------------------------------------------------------------------------------------
	285	! 6) Shift ice between categories
	286	!----------------------------------------------------------------------------------------------
	287	CALL ice_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )
	288
	289	!----------------------------------------------------------------------------------------------
	290	! 7) Make sure ht_i >= minimum ice thickness hi_min
	291	!----------------------------------------------------------------------------------------------
	292	CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1) )
	293	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) )
	294	CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) )
	295
	296	DO ji = 1, nidx
	297	IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN
	298	a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / rn_himin
	299	! MV MP 2016
	300	IF ( nn_pnd_scheme > 0 ) THEN
	301	a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin
	302	ENDIF
	303	! END MV MP 2016
	304	ht_i_1d(ji) = rn_himin
	305	ENDIF
	306	END DO
[8486]	307	!
[8422]	308	CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1) )
	309	CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) )
	310	CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) )
[8486]	311	!
[8422]	312	ENDIF
[8486]	313	!
	314	IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
	315	!
[8422]	316	END SUBROUTINE ice_itd_rem
	317
	318
	319	SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
	320	!!------------------------------------------------------------------
	321	!! * ROUTINE ice_itd_glinear *
	322	!!
	323	!! ** Purpose : build g(h) satisfying area and volume constraints (Eq. 6 and 9)
	324	!!
	325	!! ** Method : g(h) is linear and written as: g(eta) = g1(eta) + g0
	326	!! with eta = h - HL
	327	!!------------------------------------------------------------------
	328	REAL(wp), DIMENSION(:), INTENT(in ) :: HbL, HbR ! left and right category boundaries
	329	REAL(wp), DIMENSION(:), INTENT(in ) :: phice, paice ! ice thickness and concentration
	330	REAL(wp), DIMENSION(:), INTENT(inout) :: pg0, pg1 ! coefficients in linear equation for g(eta)
	331	REAL(wp), DIMENSION(:), INTENT(inout) :: phL, phR ! min and max value of range over which g(h) > 0
	332	!
	333	INTEGER :: ji ! horizontal indices
[8486]	334	REAL(wp) :: z1_3 , z2_3 ! 1/3 , 2/3
[8422]	335	REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL)
	336	REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL)
	337	REAL(wp) :: zdhr ! 1 / (hR - hL)
	338	REAL(wp) :: zwk1, zwk2 ! temporary variables
	339	!!------------------------------------------------------------------
	340	!
[8486]	341	z1_3 = 1._wp / 3._wp
	342	z2_3 = 2._wp / 3._wp
	343	!
	344	DO ji = 1, nidx
	345	!
	346	IF( paice(ji) > epsi10 .AND. phice(ji) > 0._wp ) THEN
[8422]	347	!
[8486]	348	! Initialize hL and hR
	349	phL(ji) = HbL(ji)
	350	phR(ji) = HbR(ji)
	351	!
	352	! Change hL or hR if hice falls outside central third of range,
	353	! so that hice is in the central third of the range [HL HR]
	354	zh13 = z1_3 * ( 2._wp * phL(ji) + phR(ji) )
	355	zh23 = z1_3 * ( phL(ji) + 2._wp * phR(ji) )
	356	!
	357	IF ( phice(ji) < zh13 ) THEN ; phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
	358	ELSEIF( phice(ji) > zh23 ) THEN ; phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
[8422]	359	ENDIF
	360	!
[8486]	361	! Compute coefficients of g(eta) = g0 + g1*eta
	362	zdhr = 1._wp / (phR(ji) - phL(ji))
	363	zwk1 = 6._wp * paice(ji) * zdhr
	364	zwk2 = ( phice(ji) - phL(ji) ) * zdhr
	365	pg0(ji) = zwk1 * ( z2_3 - zwk2 ) ! Eq. 14
	366	pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp ) ! Eq. 14
	367	!
	368	ELSE ! remap_flag = .false. or a_i < epsi10
	369	phL(ji) = 0._wp
	370	phR(ji) = 0._wp
	371	pg0(ji) = 0._wp
	372	pg1(ji) = 0._wp
	373	ENDIF
	374	!
	375	END DO
[8422]	376	!
	377	END SUBROUTINE ice_itd_glinear
	378
	379
	380	SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice )
	381	!!------------------------------------------------------------------
	382	!! * ROUTINE ice_itd_shiftice *
	383	!!
	384	!! ** Purpose : shift ice across category boundaries, conserving everything
	385	!! ( area, volume, energy, age*vol, and mass of salt )
	386	!!------------------------------------------------------------------
	387	INTEGER , DIMENSION(:,:), INTENT(in) :: kdonor ! donor category index
	388	REAL(wp), DIMENSION(:,:), INTENT(in) :: pdaice ! ice area transferred across boundary
	389	REAL(wp), DIMENSION(:,:), INTENT(in) :: pdvice ! ice volume transferred across boundary
[8486]	390	!
	391	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	392	INTEGER :: ii, ij, jl2, jl1 ! local integers
[8422]	393	REAL(wp) :: ztrans ! ice/snow transferred
[8486]	394	REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace
	395	REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - -
[8422]	396	!!------------------------------------------------------------------
	397
	398	CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i )
	399	CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )
	400	CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i )
	401	CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s )
	402	CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i )
	403	CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i )
	404	CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip )
	405	CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip )
	406	CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su )
	407
	408	!----------------------------------------------------------------------------------------------
	409	! 1) Define a variable equal to a_i*T_su
	410	!----------------------------------------------------------------------------------------------
	411	DO jl = 1, jpl
	412	DO ji = 1, nidx
	413	zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
	414	END DO
	415	END DO
	416
	417	!-------------------------------------------------------------------------------
	418	! 2) Transfer volume and energy between categories
	419	!-------------------------------------------------------------------------------
	420	DO jl = 1, jpl - 1
	421	DO ji = 1, nidx
[8486]	422	!
[8422]	423	jl1 = kdonor(ji,jl)
[8486]	424	IF ( jl1 == jl ) THEN ; jl2 = jl1+1
	425	ELSE ; jl2 = jl
[8422]	426	ENDIF
[8486]	427	!
	428	!!gm use of rswitch not faster as there is already IF in the DO-loop ==>>> use IF which is more readable
	429	! rswitch = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) )
	430	! zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch
	431	! !
	432	! rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) )
	433	! zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch
[8422]	434
[8486]	435	IF( v_i_2d(ji,jl1) >= epsi10 ) THEN ; zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
	436	ELSE ; zworkv(ji) = 0._wp
	437	ENDIF
	438	IF( a_i_2d(ji,jl1) >= epsi10 ) THEN ; zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
	439	ELSE ; zworka(ji) = 0._wp
	440	ENDIF
	441	!!gm end
	442	!
	443	a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas
[8422]	444	a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
[8486]	445	!
	446	v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) ! Ice volumes
[8422]	447	v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
[8486]	448	!
	449	ztrans = v_s_2d(ji,jl1) * zworkv(ji) ! Snow volumes
[8422]	450	v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
	451	v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans
[8486]	452	!
	453	! ! Ice age
	454	ztrans = oa_i_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work ????
[8422]	455	oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
	456	oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
[8486]	457	!
	458	ztrans = smv_i_2d(ji,jl1) * zworkv(ji) ! Ice salinity
	459
[8422]	460	smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - ztrans
	461	smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ztrans
[8486]	462	!
	463	! ! Surface temperature
	464	ztrans = t_su_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work ????
[8422]	465	zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans
	466	zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans
[8486]	467	!
[8422]	468	! MV MP 2016
	469	IF ( nn_pnd_scheme > 0 ) THEN
[8486]	470	! ! Pond fraction
[8422]	471	ztrans = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
	472	a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
	473	a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
[8486]	474	! ! Pond volume (also proportional to da/a)
[8422]	475	ztrans = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
	476	v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
	477	v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
	478	ENDIF
	479	! END MV MP 2016
	480	END DO
[8486]	481	!
	482	DO jk = 1, nlay_s !--- Snow heat content
	483	!
[8422]	484	DO ji = 1, nidx
	485	ii = MOD( idxice(ji) - 1, jpi ) + 1
	486	ij = ( idxice(ji) - 1 ) / jpi + 1
[8486]	487	!
[8422]	488	jl1 = kdonor(ji,jl)
	489	IF(jl1 == jl) THEN ; jl2 = jl+1
	490	ELSE ; jl2 = jl
	491	ENDIF
[8486]	492	!
[8422]	493	ztrans = e_s(ii,ij,jk,jl1) * zworkv(ji)
	494	e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
	495	e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
	496	END DO
	497	END DO
	498
[8486]	499	DO jk = 1, nlay_i !--- Ice heat content
[8422]	500	DO ji = 1, nidx
	501	ii = MOD( idxice(ji) - 1, jpi ) + 1
	502	ij = ( idxice(ji) - 1 ) / jpi + 1
[8486]	503	!
[8422]	504	jl1 = kdonor(ji,jl)
	505	IF(jl1 == jl) THEN ; jl2 = jl+1
	506	ELSE ; jl2 = jl
	507	ENDIF
[8486]	508	!
[8422]	509	ztrans = e_i(ii,ij,jk,jl1) * zworkv(ji)
	510	e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
	511	e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
	512	END DO
	513	END DO
[8486]	514	!
[8422]	515	END DO ! boundaries, 1 to jpl-1
	516
	517	!-------------------------------------------------------------------------------
	518	! 3) Update ice thickness and temperature
	519	!-------------------------------------------------------------------------------
	520	DO jl = 1, jpl
	521	DO ji = 1, nidx
	522	IF ( a_i_2d(ji,jl) > epsi10 ) THEN
[8486]	523	ht_i_2d(ji,jl) = v_i_2d(ji,jl) / a_i_2d(ji,jl)
[8422]	524	t_su_2d(ji,jl) = zaTsfn(ji,jl) / a_i_2d(ji,jl)
	525	ELSE
	526	ht_i_2d(ji,jl) = 0._wp
	527	t_su_2d(ji,jl) = rt0
	528	ENDIF
	529	END DO
	530	END DO
	531	!
[8486]	532	CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i )
	533	CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )
	534	CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i )
	535	CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s )
	536	CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i )
	537	CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i )
	538	CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip )
	539	CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip )
	540	CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su )
	541	!
[8422]	542	END SUBROUTINE ice_itd_shiftice
	543
	544
	545	SUBROUTINE ice_itd_reb
	546	!!------------------------------------------------------------------
	547	!! * ROUTINE ice_itd_reb *
	548	!!
	549	!! ** Purpose : rebin - rebins thicknesses into defined categories
	550	!!
	551	!! ** Method : If a category thickness is out of bounds, shift part (for down to top)
	552	!! or entire (for top to down) area, volume, and energy
	553	!! to the neighboring category
	554	!!------------------------------------------------------------------
	555	INTEGER :: ji, jj, jl ! dummy loop indices
	556	!
	557	INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index
	558	REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred
	559	!!------------------------------------------------------------------
[8486]	560	!
	561	jdonor(:,:) = 0
	562	zdaice(:,:) = 0._wp
	563	zdvice(:,:) = 0._wp
	564	!
	565	! !---------------------------------------
	566	DO jl = 1, jpl-1 ! identify thicknesses that are too big
	567	! !---------------------------------------
	568	nidx = 0 ; idxice(:) = 0
[8422]	569	DO jj = 1, jpj
	570	DO ji = 1, jpi
	571	IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
	572	nidx = nidx + 1
	573	idxice( nidx ) = (jj - 1) * jpi + ji
	574	ENDIF
[8486]	575	END DO
	576	END DO
	577	!
[8422]	578	!!clem CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) )
	579	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )
	580	CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )
[8486]	581	!
[8422]	582	DO ji = 1, nidx
	583	jdonor(ji,jl) = jl
	584	! how much of a_i you send in cat sup is somewhat arbitrary
	585	!!clem: these do not work properly after a restart (I do not know why)
	586	!! zdaice(ji,jl) = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji)
	587	!! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
	588	!!clem: these do not work properly after a restart (I do not know why)
	589	!! zdaice(ji,jl) = a_i_1d(ji)
	590	!! zdvice(ji,jl) = v_i_1d(ji)
	591	!!clem: these are from UCL and work ok
	592	zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp
	593	zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
	594	END DO
[8486]	595	!
[8422]	596	IF( nidx > 0 ) THEN
	597	CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl=>jl+1
	598	! Reset shift parameters
	599	jdonor(1:nidx,jl) = 0
	600	zdaice(1:nidx,jl) = 0._wp
	601	zdvice(1:nidx,jl) = 0._wp
	602	ENDIF
	603	!
	604	END DO
	605
[8486]	606	! !-----------------------------------------
	607	DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small
	608	! !-----------------------------------------
[8422]	609	nidx = 0 ; idxice(:) = 0
	610	DO jj = 1, jpj
	611	DO ji = 1, jpi
	612	IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
	613	nidx = nidx + 1
	614	idxice( nidx ) = (jj - 1) * jpi + ji
	615	ENDIF
[8486]	616	END DO
	617	END DO
	618	!
[8422]	619	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl+1) ) ! jl+1 is ok
	620	CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl+1) ) ! jl+1 is ok
	621	DO ji = 1, nidx
	622	jdonor(ji,jl) = jl + 1
	623	zdaice(ji,jl) = a_i_1d(ji)
	624	zdvice(ji,jl) = v_i_1d(ji)
	625	END DO
[8486]	626	!
[8422]	627	IF( nidx > 0 ) THEN
	628	CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl+1=>jl
	629	! Reset shift parameters
	630	jdonor(1:nidx,jl) = 0
	631	zdaice(1:nidx,jl) = 0._wp
	632	zdvice(1:nidx,jl) = 0._wp
	633	ENDIF
[8486]	634	!
[8422]	635	END DO
	636	!
	637	END SUBROUTINE ice_itd_reb
	638
[8486]	639	#else
	640	!!----------------------------------------------------------------------
	641	!! Default option : Empty module NO LIM sea-ice model
	642	!!----------------------------------------------------------------------
[8422]	643	#endif
[8486]	644
[8422]	645	!!======================================================================
	646	END MODULE iceitd

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

Download in other formats: