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

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

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

Download in other formats: