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

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

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

Download in other formats: