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 @ 8498

Last change on this file since 8498 was 8498, checked in by clem, 7 years ago
changes in style - part2 -
File size: 30.7 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	IF( v_i_2d(ji,jl1) >= epsi10 ) THEN ; zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
	429	ELSE ; zworkv(ji) = 0._wp
	430	ENDIF
	431	IF( a_i_2d(ji,jl1) >= epsi10 ) THEN ; zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
	432	ELSE ; zworka(ji) = 0._wp
	433	ENDIF
	434	!
	435	a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas
[8422]	436	a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
[8486]	437	!
	438	v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) ! Ice volumes
[8422]	439	v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
[8486]	440	!
	441	ztrans = v_s_2d(ji,jl1) * zworkv(ji) ! Snow volumes
[8422]	442	v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
	443	v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans
[8486]	444	!
	445	! ! Ice age
	446	ztrans = oa_i_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work ????
[8422]	447	oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
	448	oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
[8486]	449	!
	450	ztrans = smv_i_2d(ji,jl1) * zworkv(ji) ! Ice salinity
	451
[8422]	452	smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - ztrans
	453	smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ztrans
[8486]	454	!
	455	! ! Surface temperature
	456	ztrans = t_su_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work ????
[8422]	457	zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans
	458	zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans
[8486]	459	!
[8422]	460	! MV MP 2016
	461	IF ( nn_pnd_scheme > 0 ) THEN
[8486]	462	! ! Pond fraction
[8422]	463	ztrans = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
	464	a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
	465	a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
[8486]	466	! ! Pond volume (also proportional to da/a)
[8422]	467	ztrans = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
	468	v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
	469	v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
	470	ENDIF
	471	! END MV MP 2016
	472	END DO
[8486]	473	!
	474	DO jk = 1, nlay_s !--- Snow heat content
	475	!
[8422]	476	DO ji = 1, nidx
	477	ii = MOD( idxice(ji) - 1, jpi ) + 1
	478	ij = ( idxice(ji) - 1 ) / jpi + 1
[8486]	479	!
[8422]	480	jl1 = kdonor(ji,jl)
	481	IF(jl1 == jl) THEN ; jl2 = jl+1
	482	ELSE ; jl2 = jl
	483	ENDIF
[8486]	484	!
[8422]	485	ztrans = e_s(ii,ij,jk,jl1) * zworkv(ji)
	486	e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
	487	e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
	488	END DO
	489	END DO
	490
[8486]	491	DO jk = 1, nlay_i !--- Ice heat content
[8422]	492	DO ji = 1, nidx
	493	ii = MOD( idxice(ji) - 1, jpi ) + 1
	494	ij = ( idxice(ji) - 1 ) / jpi + 1
[8486]	495	!
[8422]	496	jl1 = kdonor(ji,jl)
	497	IF(jl1 == jl) THEN ; jl2 = jl+1
	498	ELSE ; jl2 = jl
	499	ENDIF
[8486]	500	!
[8422]	501	ztrans = e_i(ii,ij,jk,jl1) * zworkv(ji)
	502	e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
	503	e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
	504	END DO
	505	END DO
[8486]	506	!
[8422]	507	END DO ! boundaries, 1 to jpl-1
	508
	509	!-------------------------------------------------------------------------------
	510	! 3) Update ice thickness and temperature
	511	!-------------------------------------------------------------------------------
	512	DO jl = 1, jpl
	513	DO ji = 1, nidx
	514	IF ( a_i_2d(ji,jl) > epsi10 ) THEN
[8486]	515	ht_i_2d(ji,jl) = v_i_2d(ji,jl) / a_i_2d(ji,jl)
[8422]	516	t_su_2d(ji,jl) = zaTsfn(ji,jl) / a_i_2d(ji,jl)
	517	ELSE
	518	ht_i_2d(ji,jl) = 0._wp
	519	t_su_2d(ji,jl) = rt0
	520	ENDIF
	521	END DO
	522	END DO
	523	!
[8486]	524	CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i )
	525	CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i )
	526	CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i )
	527	CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s )
	528	CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i )
	529	CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i )
	530	CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip )
	531	CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip )
	532	CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su )
	533	!
[8422]	534	END SUBROUTINE ice_itd_shiftice
	535
	536
	537	SUBROUTINE ice_itd_reb
	538	!!------------------------------------------------------------------
	539	!! * ROUTINE ice_itd_reb *
	540	!!
	541	!! ** Purpose : rebin - rebins thicknesses into defined categories
	542	!!
	543	!! ** Method : If a category thickness is out of bounds, shift part (for down to top)
	544	!! or entire (for top to down) area, volume, and energy
	545	!! to the neighboring category
	546	!!------------------------------------------------------------------
	547	INTEGER :: ji, jj, jl ! dummy loop indices
	548	!
	549	INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index
	550	REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred
	551	!!------------------------------------------------------------------
[8486]	552	!
	553	jdonor(:,:) = 0
	554	zdaice(:,:) = 0._wp
	555	zdvice(:,:) = 0._wp
	556	!
	557	! !---------------------------------------
	558	DO jl = 1, jpl-1 ! identify thicknesses that are too big
	559	! !---------------------------------------
	560	nidx = 0 ; idxice(:) = 0
[8422]	561	DO jj = 1, jpj
	562	DO ji = 1, jpi
	563	IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
	564	nidx = nidx + 1
	565	idxice( nidx ) = (jj - 1) * jpi + ji
	566	ENDIF
[8486]	567	END DO
	568	END DO
	569	!
[8422]	570	!!clem CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) )
	571	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) )
	572	CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) )
[8486]	573	!
[8422]	574	DO ji = 1, nidx
	575	jdonor(ji,jl) = jl
	576	! how much of a_i you send in cat sup is somewhat arbitrary
	577	!!clem: these do not work properly after a restart (I do not know why)
	578	!! zdaice(ji,jl) = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji)
	579	!! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
	580	!!clem: these do not work properly after a restart (I do not know why)
	581	!! zdaice(ji,jl) = a_i_1d(ji)
	582	!! zdvice(ji,jl) = v_i_1d(ji)
	583	!!clem: these are from UCL and work ok
	584	zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp
	585	zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
	586	END DO
[8486]	587	!
[8422]	588	IF( nidx > 0 ) THEN
	589	CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl=>jl+1
	590	! Reset shift parameters
	591	jdonor(1:nidx,jl) = 0
	592	zdaice(1:nidx,jl) = 0._wp
	593	zdvice(1:nidx,jl) = 0._wp
	594	ENDIF
	595	!
	596	END DO
	597
[8486]	598	! !-----------------------------------------
	599	DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small
	600	! !-----------------------------------------
[8422]	601	nidx = 0 ; idxice(:) = 0
	602	DO jj = 1, jpj
	603	DO ji = 1, jpi
	604	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
	605	nidx = nidx + 1
	606	idxice( nidx ) = (jj - 1) * jpi + ji
	607	ENDIF
[8486]	608	END DO
	609	END DO
	610	!
[8422]	611	CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl+1) ) ! jl+1 is ok
	612	CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl+1) ) ! jl+1 is ok
	613	DO ji = 1, nidx
	614	jdonor(ji,jl) = jl + 1
	615	zdaice(ji,jl) = a_i_1d(ji)
	616	zdvice(ji,jl) = v_i_1d(ji)
	617	END DO
[8486]	618	!
[8422]	619	IF( nidx > 0 ) THEN
	620	CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl+1=>jl
	621	! Reset shift parameters
	622	jdonor(1:nidx,jl) = 0
	623	zdaice(1:nidx,jl) = 0._wp
	624	zdvice(1:nidx,jl) = 0._wp
	625	ENDIF
[8486]	626	!
[8422]	627	END DO
	628	!
	629	END SUBROUTINE ice_itd_reb
	630
[8486]	631	#else
	632	!!----------------------------------------------------------------------
	633	!! Default option : Empty module NO LIM sea-ice model
	634	!!----------------------------------------------------------------------
[8422]	635	#endif
[8486]	636
[8422]	637	!!======================================================================
	638	END MODULE iceitd

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

Download in other formats: