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 NEMO/releases/r4.0/r4.0-HEAD/src/ICE – NEMO

Context Navigation

source: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/iceitd.F90 @ 14675

Last change on this file since 14675 was 14026, checked in by clem, 4 years ago
4.0-HEAD: fix bugs and defects related to tickets #2573 #2575 #2576 #2578. Sette passed and those fixes are now in the trunk, so unless there is a tricky trick somewhere, everything should be fine.
Property svn:keywords set to `Id`
File size: 37.9 KB

Rev	Line
[8586]	1	MODULE iceitd
	2	!!======================================================================
	3	!! * MODULE iceitd *
	4	!! sea-ice : ice thickness distribution
	5	!!======================================================================
[9604]	6	!! History : 3.0 ! 2005-12 (M. Vancoppenolle) original code (based on CICE)
	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	8	!!----------------------------------------------------------------------
[9570]	9	#if defined key_si3
[8586]	10	!!----------------------------------------------------------------------
[9570]	11	!! 'key_si3' SI3 sea-ice model
[8586]	12	!!----------------------------------------------------------------------
[8813]	13	!! ice_itd_rem : redistribute ice thicknesses after thermo growth and melt
	14	!! itd_glinear : build g(h) satisfying area and volume constraints
	15	!! itd_shiftice : shift ice across category boundaries, conserving everything
	16	!! ice_itd_reb : rebin ice thicknesses into bounded categories
	17	!! ice_itd_init : read ice thicknesses mean and min from namelist
[8586]	18	!!----------------------------------------------------------------------
	19	USE dom_oce ! ocean domain
	20	USE phycst ! physical constants
	21	USE ice1D ! sea-ice: thermodynamic variables
	22	USE ice ! sea-ice: variables
[10994]	23	USE icevar ! sea-ice: operations
[8586]	24	USE icectl ! sea-ice: conservation tests
	25	USE icetab ! sea-ice: convert 1D<=>2D
	26	!
	27	USE in_out_manager ! I/O manager
	28	USE lib_mpp ! MPP library
	29	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	30	USE prtctl ! Print control
[14026]	31	USE timing ! Timing
[8586]	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC ice_itd_init ! called in icestp
	37	PUBLIC ice_itd_rem ! called in icethd
	38	PUBLIC ice_itd_reb ! called in icecor
	39
[8813]	40	INTEGER :: nice_catbnd ! choice of the type of ice category function
	41	! ! associated indices:
	42	INTEGER, PARAMETER :: np_cathfn = 1 ! categories defined by a function
	43	INTEGER, PARAMETER :: np_catusr = 2 ! categories defined by the user
	44	!
	45	! !! namelist (namitd)
	46	LOGICAL :: ln_cat_hfn ! ice categories are defined by function like rn_himean**(-0.05)
	47	REAL(wp) :: rn_himean ! mean thickness of the domain
	48	LOGICAL :: ln_cat_usr ! ice categories are defined by rn_catbnd
	49	REAL(wp), DIMENSION(0:100) :: rn_catbnd ! ice categories bounds
[13284]	50	REAL(wp) :: rn_himax ! maximum ice thickness allowed
[8813]	51	!
[8586]	52	!!----------------------------------------------------------------------
[9598]	53	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	54	!! $Id$
[10068]	55	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	56	!!----------------------------------------------------------------------
	57	CONTAINS
	58
	59	SUBROUTINE ice_itd_rem( kt )
	60	!!------------------------------------------------------------------
	61	!! * ROUTINE ice_itd_rem *
	62	!!
	63	!! ** Purpose : computes the redistribution of ice thickness
	64	!! after thermodynamic growth of ice thickness
	65	!!
	66	!! ** Method : Linear remapping
	67	!!
	68	!! References : W.H. Lipscomb, JGR 2001
	69	!!------------------------------------------------------------------
	70	INTEGER , INTENT (in) :: kt ! Ocean time step
	71	!
	72	INTEGER :: ji, jj, jl, jcat ! dummy loop index
	73	INTEGER :: ipti ! local integer
	74	REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars
	75	REAL(wp) :: zx2, zwk2, zda0, zetamax ! - -
	76	REAL(wp) :: zx3
	77	REAL(wp) :: zslope ! used to compute local thermodynamic "speeds"
[8813]	78	!
[8586]	79	INTEGER , DIMENSION(jpij) :: iptidx ! compute remapping or not
	80	INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index
	81	REAL(wp), DIMENSION(jpij,jpl) :: zdhice ! ice thickness increment
	82	REAL(wp), DIMENSION(jpij,jpl) :: g0, g1 ! coefficients for fitting the line of the ITD
	83	REAL(wp), DIMENSION(jpij,jpl) :: hL, hR ! left and right boundary for the ITD for each thickness
	84	REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! local increment of ice area and volume
	85	REAL(wp), DIMENSION(jpij) :: zhb0, zhb1 ! category boundaries for thinnes categories
	86	REAL(wp), DIMENSION(jpij,0:jpl) :: zhbnew ! new boundaries of ice categories
	87	!!------------------------------------------------------------------
[14026]	88	IF( ln_timing ) CALL timing_start('iceitd_rem')
[8586]	89
	90	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution'
	91
	92	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[11536]	93	IF( ln_icediachk ) CALL ice_cons2D (0, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
[8586]	94
	95	!-----------------------------------------------------------------------------------------------
	96	! 1) Identify grid cells with ice
	97	!-----------------------------------------------------------------------------------------------
[10994]	98	at_i(:,:) = SUM( a_i, dim=3 )
	99	!
[8586]	100	npti = 0 ; nptidx(:) = 0
	101	DO jj = 1, jpj
	102	DO ji = 1, jpi
	103	IF ( at_i(ji,jj) > epsi10 ) THEN
[9880]	104	npti = npti + 1
[8586]	105	nptidx( npti ) = (jj - 1) * jpi + ji
	106	ENDIF
	107	END DO
	108	END DO
	109
	110	!-----------------------------------------------------------------------------------------------
	111	! 2) Compute new category boundaries
	112	!-----------------------------------------------------------------------------------------------
	113	IF( npti > 0 ) THEN
[8813]	114	!
[8586]	115	zdhice(:,:) = 0._wp
	116	zhbnew(:,:) = 0._wp
[8813]	117	!
[8586]	118	CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i )
	119	CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
[9880]	120	CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i )
	121	CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d(1:npti,1:jpl), a_i_b )
[8813]	122	!
[8586]	123	DO jl = 1, jpl
	124	! Compute thickness change in each ice category
	125	DO ji = 1, npti
[9880]	126	IF( a_i_2d(ji,jl) > epsi10 ) zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)
[8586]	127	END DO
	128	END DO
[8813]	129	!
[8586]	130	! --- New boundaries for category 1:jpl-1 --- !
	131	DO jl = 1, jpl - 1
	132	!
	133	DO ji = 1, npti
	134	!
	135	! --- New boundary: Hn* = Hn + Fn*dt --- !
	136	! Fndt = ( fn + (fn+1 - fn)/(hn+1 - hn) (Hn - hn) ) * dt = zdhice + zslope * (Hmax - h_i_b)
	137	!
	138	IF ( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl+1) & a(jl) /= 0
[9880]	139	zslope = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( h_ib_2d(ji,jl+1) - h_ib_2d(ji,jl) )
[8586]	140	zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - h_ib_2d(ji,jl) )
	141	ELSEIF( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN ! a(jl+1)=0 => Hn* = Hn + fn*dt
	142	zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl)
	143	ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl)=0 => Hn* = Hn + fn+1*dt
	144	zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1)
	145	ELSE ! a(jl+1) & a(jl) = 0
	146	zhbnew(ji,jl) = hi_max(jl)
	147	ENDIF
	148	!
	149	! --- 2 conditions for remapping --- !
	150	! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi
	151	! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
[8813]	152	! in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
[8586]	153	IF( a_i_2d(ji,jl ) > epsi10 .AND. h_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) nptidx(ji) = 0
	154	IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) nptidx(ji) = 0
[8813]	155	!
[8586]	156	! 2) Hn-1 < Hn* < Hn+1
	157	IF( zhbnew(ji,jl) < hi_max(jl-1) ) nptidx(ji) = 0
	158	IF( zhbnew(ji,jl) > hi_max(jl+1) ) nptidx(ji) = 0
[8813]	159	!
[8586]	160	END DO
	161	END DO
	162	!
	163	! --- New boundaries for category jpl --- !
	164	DO ji = 1, npti
	165	IF( a_i_2d(ji,jpl) > epsi10 ) THEN
	166	zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
	167	ELSE
	168	zhbnew(ji,jpl) = hi_max(jpl)
	169	ENDIF
[8813]	170	!
[8586]	171	! --- 1 additional condition for remapping (1st category) --- !
	172	! H0+epsi < h1(t) < H1-epsi
	173	! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible
[8813]	174	! in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
[8586]	175	IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) nptidx(ji) = 0
	176	IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) nptidx(ji) = 0
	177	END DO
	178	!
	179	!-----------------------------------------------------------------------------------------------
	180	! 3) Identify cells where remapping
	181	!-----------------------------------------------------------------------------------------------
[8813]	182	ipti = 0 ; iptidx(:) = 0
[8586]	183	DO ji = 1, npti
	184	IF( nptidx(ji) /= 0 ) THEN
	185	ipti = ipti + 1
	186	iptidx(ipti) = nptidx(ji)
	187	zhbnew(ipti,:) = zhbnew(ji,:) ! adjust zhbnew to new indices
	188	ENDIF
	189	END DO
	190	nptidx(:) = iptidx(:)
	191	npti = ipti
	192	!
	193	ENDIF
	194
	195	!-----------------------------------------------------------------------------------------------
	196	! 4) Compute g(h)
	197	!-----------------------------------------------------------------------------------------------
	198	IF( npti > 0 ) THEN
	199	!
	200	zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1)
	201	g0(:,:) = 0._wp ; g1(:,:) = 0._wp
	202	hL(:,:) = 0._wp ; hR(:,:) = 0._wp
	203	!
	204	DO jl = 1, jpl
	205	!
	206	CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) )
[9880]	207	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) )
	208	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) )
	209	CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,jl) )
[8586]	210	!
	211	IF( jl == 1 ) THEN
	212	!
	213	! --- g(h) for category 1 --- !
[8813]	214	CALL itd_glinear( zhb0(1:npti) , zhb1(1:npti) , h_ib_1d(1:npti) , a_i_1d(1:npti) , & ! in
[9880]	215	& g0 (1:npti,1), g1 (1:npti,1), hL (1:npti,1), hR (1:npti,1) ) ! out
[11732]	216	!
[8586]	217	! Area lost due to melting of thin ice
	218	DO ji = 1, npti
	219	!
	220	IF( a_i_1d(ji) > epsi10 ) THEN
	221	!
	222	zdh0 = h_i_1d(ji) - h_ib_1d(ji)
[11732]	223	IF( zdh0 < 0.0 ) THEN ! remove area from category 1
[8586]	224	zdh0 = MIN( -zdh0, hi_max(1) )
	225	!Integrate g(1) from 0 to dh0 to estimate area melted
	226	zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1)
	227	!
	228	IF( zetamax > 0.0 ) THEN
	229	zx1 = zetamax
	230	zx2 = 0.5 * zetamax * zetamax
[11732]	231	zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed
[8586]	232	zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i
[11732]	233	zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting of thin ice (zdamax > 0)
[8586]	234	! Remove area, conserving volume
	235	h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
[9880]	236	a_i_1d(ji) = a_i_1d(ji) - zda0
	237	v_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) ! useless ?
[8586]	238	ENDIF
	239	!
	240	ELSE ! if ice accretion zdh0 > 0
	241	! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
	242	zhbnew(ji,0) = MIN( zdh0, hi_max(1) )
	243	ENDIF
	244	!
	245	ENDIF
	246	!
	247	END DO
	248	!
[9880]	249	CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
	250	CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
	251	CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
[8586]	252	!
	253	ENDIF ! jl=1
	254	!
	255	! --- g(h) for each thickness category --- !
[8813]	256	CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti) , a_i_1d(1:npti) , & ! in
[10994]	257	& g0 (1:npti,jl ), g1 (1:npti,jl), hL (1:npti,jl), hR (1:npti,jl) ) ! out
[8586]	258	!
	259	END DO
	260
	261	!-----------------------------------------------------------------------------------------------
	262	! 5) Compute area and volume to be shifted across each boundary (Eq. 18)
	263	!-----------------------------------------------------------------------------------------------
	264	DO jl = 1, jpl - 1
	265	!
	266	DO ji = 1, npti
	267	!
	268	! left and right integration limits in eta space
	269	IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1
	270	zetamin = MAX( hi_max(jl) , hL(ji,jl) ) - hL(ji,jl) ! hi_max(jl) - hL
	271	zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl) ! hR - hL
	272	jdonor(ji,jl) = jl
	273	ELSE ! Hn* <= Hn => transfer from jl+1 to jl
	274	zetamin = 0.0
	275	zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1) ! hi_max(jl) - hL
	276	jdonor(ji,jl) = jl + 1
	277	ENDIF
	278	zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin
	279	!
	280	zx1 = zetamax - zetamin
	281	zwk1 = zetamin * zetamin
	282	zwk2 = zetamax * zetamax
	283	zx2 = 0.5 * ( zwk2 - zwk1 )
	284	zwk1 = zwk1 * zetamin
	285	zwk2 = zwk2 * zetamax
	286	zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 )
	287	jcat = jdonor(ji,jl)
	288	zdaice(ji,jl) = g1(ji,jcat)zx2 + g0(ji,jcat)zx1
	289	zdvice(ji,jl) = g1(ji,jcat)zx3 + g0(ji,jcat)zx2 + zdaice(ji,jl)*hL(ji,jcat)
	290	!
	291	END DO
	292	END DO
	293
	294	!----------------------------------------------------------------------------------------------
	295	! 6) Shift ice between categories
	296	!----------------------------------------------------------------------------------------------
[8813]	297	CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
[8586]	298
	299	!----------------------------------------------------------------------------------------------
	300	! 7) Make sure h_i >= minimum ice thickness hi_min
	301	!----------------------------------------------------------------------------------------------
[8906]	302	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
	303	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
	304	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
[8813]	305	!
[8586]	306	DO ji = 1, npti
	307	IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
[9880]	308	a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin
[13284]	309	IF( ln_pnd_LEV ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
[8586]	310	h_i_1d(ji) = rn_himin
	311	ENDIF
	312	END DO
	313	!
[8906]	314	CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
	315	CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
	316	CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
[8586]	317	!
	318	ENDIF
	319	!
	320	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[11536]	321	IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
[14026]	322	IF( ln_timing ) CALL timing_stop ('iceitd_rem')
[8586]	323	!
	324	END SUBROUTINE ice_itd_rem
	325
	326
[8813]	327	SUBROUTINE itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR )
[8586]	328	!!------------------------------------------------------------------
[8813]	329	!! * ROUTINE itd_glinear *
[8586]	330	!!
	331	!! ** Purpose : build g(h) satisfying area and volume constraints (Eq. 6 and 9)
	332	!!
	333	!! ** Method : g(h) is linear and written as: g(eta) = g1(eta) + g0
	334	!! with eta = h - HL
	335	!!------------------------------------------------------------------
	336	REAL(wp), DIMENSION(:), INTENT(in ) :: HbL, HbR ! left and right category boundaries
	337	REAL(wp), DIMENSION(:), INTENT(in ) :: phice, paice ! ice thickness and concentration
	338	REAL(wp), DIMENSION(:), INTENT(inout) :: pg0, pg1 ! coefficients in linear equation for g(eta)
	339	REAL(wp), DIMENSION(:), INTENT(inout) :: phL, phR ! min and max value of range over which g(h) > 0
	340	!
	341	INTEGER :: ji ! horizontal indices
	342	REAL(wp) :: z1_3 , z2_3 ! 1/3 , 2/3
	343	REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL)
	344	REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL)
	345	REAL(wp) :: zdhr ! 1 / (hR - hL)
	346	REAL(wp) :: zwk1, zwk2 ! temporary variables
	347	!!------------------------------------------------------------------
	348	!
	349	z1_3 = 1._wp / 3._wp
	350	z2_3 = 2._wp / 3._wp
	351	!
	352	DO ji = 1, npti
	353	!
[11732]	354	IF( paice(ji) > epsi10 .AND. phice(ji) > epsi10 ) THEN
[8586]	355	!
	356	! Initialize hL and hR
	357	phL(ji) = HbL(ji)
	358	phR(ji) = HbR(ji)
	359	!
	360	! Change hL or hR if hice falls outside central third of range,
	361	! so that hice is in the central third of the range [HL HR]
	362	zh13 = z1_3 * ( 2._wp * phL(ji) + phR(ji) )
	363	zh23 = z1_3 * ( phL(ji) + 2._wp * phR(ji) )
	364	!
	365	IF ( phice(ji) < zh13 ) THEN ; phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left
	366	ELSEIF( phice(ji) > zh23 ) THEN ; phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right
	367	ENDIF
	368	!
	369	! Compute coefficients of g(eta) = g0 + g1*eta
	370	zdhr = 1._wp / (phR(ji) - phL(ji))
	371	zwk1 = 6._wp * paice(ji) * zdhr
	372	zwk2 = ( phice(ji) - phL(ji) ) * zdhr
	373	pg0(ji) = zwk1 * ( z2_3 - zwk2 ) ! Eq. 14
	374	pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5_wp ) ! Eq. 14
	375	!
	376	ELSE ! remap_flag = .false. or a_i < epsi10
	377	phL(ji) = 0._wp
	378	phR(ji) = 0._wp
	379	pg0(ji) = 0._wp
	380	pg1(ji) = 0._wp
	381	ENDIF
	382	!
	383	END DO
	384	!
[8813]	385	END SUBROUTINE itd_glinear
[8586]	386
	387
[8813]	388	SUBROUTINE itd_shiftice( kdonor, pdaice, pdvice )
[8586]	389	!!------------------------------------------------------------------
[8813]	390	!! * ROUTINE itd_shiftice *
[8586]	391	!!
	392	!! ** Purpose : shift ice across category boundaries, conserving everything
	393	!! ( area, volume, energy, age*vol, and mass of salt )
	394	!!------------------------------------------------------------------
	395	INTEGER , DIMENSION(:,:), INTENT(in) :: kdonor ! donor category index
	396	REAL(wp), DIMENSION(:,:), INTENT(in) :: pdaice ! ice area transferred across boundary
	397	REAL(wp), DIMENSION(:,:), INTENT(in) :: pdvice ! ice volume transferred across boundary
	398	!
[10994]	399	INTEGER :: ji, jl, jk ! dummy loop indices
	400	INTEGER :: jl2, jl1 ! local integers
[8586]	401	REAL(wp) :: ztrans ! ice/snow transferred
[10994]	402	REAL(wp), DIMENSION(jpij) :: zworka, zworkv ! workspace
	403	REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn ! - -
	404	REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d
	405	REAL(wp), DIMENSION(jpij,nlay_s,jpl) :: ze_s_2d
[8586]	406	!!------------------------------------------------------------------
	407
	408	CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i )
	409	CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i )
	410	CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i )
	411	CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s )
	412	CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
	413	CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
	414	CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
	415	CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
[13284]	416	CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
[8586]	417	CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
[10994]	418	DO jl = 1, jpl
	419	DO jk = 1, nlay_s
	420	CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
	421	END DO
	422	DO jk = 1, nlay_i
	423	CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
	424	END DO
	425	END DO
	426	! to correct roundoff errors on a_i
	427	CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti), rn_amax_2d )
[8586]	428
	429	!----------------------------------------------------------------------------------------------
	430	! 1) Define a variable equal to a_i*T_su
	431	!----------------------------------------------------------------------------------------------
	432	DO jl = 1, jpl
	433	DO ji = 1, npti
	434	zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
	435	END DO
	436	END DO
	437
	438	!-------------------------------------------------------------------------------
	439	! 2) Transfer volume and energy between categories
	440	!-------------------------------------------------------------------------------
	441	DO jl = 1, jpl - 1
	442	DO ji = 1, npti
	443	!
	444	jl1 = kdonor(ji,jl)
	445	!
	446	IF( jl1 > 0 ) THEN
	447	!
	448	IF ( jl1 == jl ) THEN ; jl2 = jl1+1
	449	ELSE ; jl2 = jl
	450	ENDIF
	451	!
	452	IF( v_i_2d(ji,jl1) >= epsi10 ) THEN ; zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
	453	ELSE ; zworkv(ji) = 0._wp
	454	ENDIF
	455	IF( a_i_2d(ji,jl1) >= epsi10 ) THEN ; zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
	456	ELSE ; zworka(ji) = 0._wp
	457	ENDIF
	458	!
	459	a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) ! Ice areas
	460	a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
	461	!
	462	v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) ! Ice volumes
	463	v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
	464	!
	465	ztrans = v_s_2d(ji,jl1) * zworkv(ji) ! Snow volumes
	466	v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
	467	v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans
[9880]	468	!
	469	ztrans = oa_i_2d(ji,jl1) * zworka(ji) ! Ice age
[8586]	470	oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans
	471	oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans
	472	!
	473	ztrans = sv_i_2d(ji,jl1) * zworkv(ji) ! Ice salinity
	474	sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - ztrans
	475	sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ztrans
	476	!
[9880]	477	ztrans = zaTsfn(ji,jl1) * zworka(ji) ! Surface temperature
[8586]	478	zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans
	479	zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans
	480	!
[13284]	481	IF ( ln_pnd_LEV ) THEN
[9880]	482	ztrans = a_ip_2d(ji,jl1) * zworka(ji) ! Pond fraction
[8586]	483	a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
	484	a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
[9880]	485	!
[14026]	486	!!$ ztrans = v_ip_2d(ji,jl1) * zworka(ji) ! Pond volume (also proportional to da/a)
	487	ztrans = v_ip_2d(ji,jl1) * zworkv(ji) ! Pond volume
[8586]	488	v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
	489	v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
[13284]	490	!
	491	IF ( ln_pnd_lids ) THEN ! Pond lid volume
[14026]	492	!!$ ztrans = v_il_2d(ji,jl1) * zworka(ji)
	493	ztrans = v_il_2d(ji,jl1) * zworkv(ji)
[13284]	494	v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans
	495	v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans
	496	ENDIF
[8586]	497	ENDIF
	498	!
	499	ENDIF ! jl1 >0
	500	END DO
	501	!
	502	DO jk = 1, nlay_s !--- Snow heat content
	503	DO ji = 1, npti
	504	!
	505	jl1 = kdonor(ji,jl)
	506	!
	507	IF( jl1 > 0 ) THEN
	508	IF(jl1 == jl) THEN ; jl2 = jl+1
	509	ELSE ; jl2 = jl
	510	ENDIF
[10994]	511	ztrans = ze_s_2d(ji,jk,jl1) * zworkv(ji)
	512	ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) - ztrans
	513	ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ztrans
[8586]	514	ENDIF
	515	END DO
	516	END DO
[8813]	517	!
[8586]	518	DO jk = 1, nlay_i !--- Ice heat content
	519	DO ji = 1, npti
	520	!
	521	jl1 = kdonor(ji,jl)
	522	!
	523	IF( jl1 > 0 ) THEN
	524	IF(jl1 == jl) THEN ; jl2 = jl+1
	525	ELSE ; jl2 = jl
	526	ENDIF
[10994]	527	ztrans = ze_i_2d(ji,jk,jl1) * zworkv(ji)
	528	ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) - ztrans
	529	ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + ztrans
[8586]	530	ENDIF
	531	END DO
	532	END DO
	533	!
	534	END DO ! boundaries, 1 to jpl-1
[10994]	535
	536	!-------------------
	537	! 3) roundoff errors
	538	!-------------------
	539	! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
	540	! because of truncation error ( i.e. 1. - 1. /= 0 )
[13284]	541	CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d )
[10994]	542
	543	! at_i must be <= rn_amax
	544	zworka(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
	545	DO jl = 1, jpl
	546	WHERE( zworka(1:npti) > rn_amax_1d(1:npti) ) &
	547	& a_i_2d(1:npti,jl) = a_i_2d(1:npti,jl) * rn_amax_1d(1:npti) / zworka(1:npti)
	548	END DO
[8586]	549
	550	!-------------------------------------------------------------------------------
[10994]	551	! 4) Update ice thickness and temperature
[8586]	552	!-------------------------------------------------------------------------------
	553	WHERE( a_i_2d(1:npti,:) >= epsi20 )
[9880]	554	h_i_2d (1:npti,:) = v_i_2d(1:npti,:) / a_i_2d(1:npti,:)
[8586]	555	t_su_2d(1:npti,:) = zaTsfn(1:npti,:) / a_i_2d(1:npti,:)
	556	ELSEWHERE
[9880]	557	h_i_2d (1:npti,:) = 0._wp
[8586]	558	t_su_2d(1:npti,:) = rt0
	559	END WHERE
	560	!
	561	CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i )
	562	CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i )
	563	CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i )
	564	CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s )
	565	CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
	566	CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
	567	CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
	568	CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
[13284]	569	CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il )
[8586]	570	CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
[10994]	571	DO jl = 1, jpl
	572	DO jk = 1, nlay_s
	573	CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
	574	END DO
	575	DO jk = 1, nlay_i
	576	CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
	577	END DO
	578	END DO
[8586]	579	!
[8813]	580	END SUBROUTINE itd_shiftice
[8586]	581
	582
	583	SUBROUTINE ice_itd_reb( kt )
	584	!!------------------------------------------------------------------
	585	!! * ROUTINE ice_itd_reb *
	586	!!
	587	!! ** Purpose : rebin - rebins thicknesses into defined categories
	588	!!
	589	!! ** Method : If a category thickness is out of bounds, shift part (for down to top)
	590	!! or entire (for top to down) area, volume, and energy
	591	!! to the neighboring category
	592	!!------------------------------------------------------------------
	593	INTEGER , INTENT (in) :: kt ! Ocean time step
	594	INTEGER :: ji, jj, jl ! dummy loop indices
	595	!
	596	INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index
	597	REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred
	598	!!------------------------------------------------------------------
[14026]	599	IF( ln_timing ) CALL timing_start('iceitd_reb')
[8586]	600	!
	601	IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'
[8813]	602	!
[10994]	603	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[11536]	604	IF( ln_icediachk ) CALL ice_cons2D (0, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
[10994]	605	!
[8586]	606	jdonor(:,:) = 0
	607	zdaice(:,:) = 0._wp
	608	zdvice(:,:) = 0._wp
	609	!
	610	! !---------------------------------------
	611	DO jl = 1, jpl-1 ! identify thicknesses that are too big
	612	! !---------------------------------------
	613	npti = 0 ; nptidx(:) = 0
	614	DO jj = 1, jpj
	615	DO ji = 1, jpi
[9880]	616	IF( a_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
[8586]	617	npti = npti + 1
	618	nptidx( npti ) = (jj - 1) * jpi + ji
	619	ENDIF
	620	END DO
	621	END DO
	622	!
[13617]	623	IF( npti > 0 ) THEN
	624	!!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
	625	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
	626	CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
	627	!
	628	DO ji = 1, npti
	629	jdonor(ji,jl) = jl
	630	! how much of a_i you send in cat sup is somewhat arbitrary
[14026]	631	! these are from CICE => transfer everything
	632	!!zdaice(ji,jl) = a_i_1d(ji)
	633	!!zdvice(ji,jl) = v_i_1d(ji)
	634	! these are from LLN => transfer only half of the category
	635	zdaice(ji,jl) = 0.5_wp * a_i_1d(ji)
	636	zdvice(ji,jl) = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl)
[13617]	637	END DO
	638	!
[8813]	639	CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl=>jl+1
[8586]	640	! Reset shift parameters
	641	jdonor(1:npti,jl) = 0
	642	zdaice(1:npti,jl) = 0._wp
	643	zdvice(1:npti,jl) = 0._wp
	644	ENDIF
	645	!
	646	END DO
	647
	648	! !-----------------------------------------
	649	DO jl = jpl-1, 1, -1 ! Identify thicknesses that are too small
	650	! !-----------------------------------------
	651	npti = 0 ; nptidx(:) = 0
	652	DO jj = 1, jpj
	653	DO ji = 1, jpi
[9880]	654	IF( a_i(ji,jj,jl+1) > 0._wp .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
[8586]	655	npti = npti + 1
	656	nptidx( npti ) = (jj - 1) * jpi + ji
	657	ENDIF
	658	END DO
	659	END DO
	660	!
	661	IF( npti > 0 ) THEN
[13617]	662	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
	663	CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
	664	!
	665	DO ji = 1, npti
	666	jdonor(ji,jl) = jl + 1
	667	zdaice(ji,jl) = a_i_1d(ji)
	668	zdvice(ji,jl) = v_i_1d(ji)
	669	END DO
	670	!
[8813]	671	CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl+1=>jl
[8586]	672	! Reset shift parameters
	673	jdonor(1:npti,jl) = 0
	674	zdaice(1:npti,jl) = 0._wp
	675	zdvice(1:npti,jl) = 0._wp
	676	ENDIF
	677	!
	678	END DO
	679	!
[10994]	680	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[11536]	681	IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
[14026]	682	IF( ln_timing ) CALL timing_stop ('iceitd_reb')
[10994]	683	!
[8586]	684	END SUBROUTINE ice_itd_reb
	685
[8813]	686
[8586]	687	SUBROUTINE ice_itd_init
	688	!!------------------------------------------------------------------
	689	!! * ROUTINE ice_itd_init *
	690	!!
	691	!! ** Purpose : Initializes the ice thickness distribution
	692	!! ** Method : ...
	693	!! ** input : Namelist namitd
	694	!!-------------------------------------------------------------------
[8813]	695	INTEGER :: jl ! dummy loop index
	696	INTEGER :: ios, ioptio ! Local integer output status for namelist read
[8586]	697	REAL(wp) :: zhmax, znum, zden, zalpha ! - -
[8813]	698	!
[13284]	699	NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin, rn_himax
[8586]	700	!!------------------------------------------------------------------
	701	!
	702	REWIND( numnam_ice_ref ) ! Namelist namitd in reference namelist : Parameters for ice
	703	READ ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
[11536]	704	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist' )
[8586]	705	REWIND( numnam_ice_cfg ) ! Namelist namitd in configuration namelist : Parameters for ice
	706	READ ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
[11536]	707	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namitd in configuration namelist' )
[9169]	708	IF(lwm) WRITE( numoni, namitd )
[8586]	709	!
	710	IF(lwp) THEN ! control print
	711	WRITE(numout,*)
	712	WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
	713	WRITE(numout,*) '~~~~~~~~~~~~'
	714	WRITE(numout,*) ' Namelist namitd: '
[8813]	715	WRITE(numout,) ' Ice categories are defined by a function of rn_himean*(-0.05) ln_cat_hfn = ', ln_cat_hfn
	716	WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean
	717	WRITE(numout,*) ' Ice categories are defined by rn_catbnd ln_cat_usr = ', ln_cat_usr
[13284]	718	WRITE(numout,*) ' minimum ice thickness allowed rn_himin = ', rn_himin
	719	WRITE(numout,*) ' maximum ice thickness allowed rn_himax = ', rn_himax
[8586]	720	ENDIF
	721	!
	722	!-----------------------------------!
	723	! Thickness categories boundaries !
	724	!-----------------------------------!
[8813]	725	! !== set the choice of ice categories ==!
	726	ioptio = 0
	727	IF( ln_cat_hfn ) THEN ; ioptio = ioptio + 1 ; nice_catbnd = np_cathfn ; ENDIF
	728	IF( ln_cat_usr ) THEN ; ioptio = ioptio + 1 ; nice_catbnd = np_catusr ; ENDIF
	729	IF( ioptio /= 1 ) CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' )
[8586]	730	!
[8813]	731	SELECT CASE( nice_catbnd )
	732	! !------------------------!
	733	CASE( np_cathfn ) ! h^(-alpha) function
	734	! !------------------------!
	735	zalpha = 0.05_wp
	736	zhmax = 3._wp * rn_himean
[8966]	737	hi_max(0) = 0._wp
[8813]	738	DO jl = 1, jpl
	739	znum = jpl * ( zhmax+1 )**zalpha
	740	zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
	741	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
	742	END DO
	743	! !------------------------!
	744	CASE( np_catusr ) ! user defined
	745	! !------------------------!
	746	DO jl = 0, jpl
	747	hi_max(jl) = rn_catbnd(jl)
	748	END DO
	749	!
	750	END SELECT
[8586]	751	!
	752	DO jl = 1, jpl ! mean thickness by category
	753	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
	754	END DO
	755	!
[13284]	756	hi_max(jpl) = rn_himax ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
[8586]	757	!
	758	IF(lwp) WRITE(numout,*)
	759	IF(lwp) WRITE(numout,*) ' ===>>> resulting thickness category boundaries :'
	760	IF(lwp) WRITE(numout,*) ' hi_max(:)= ', hi_max(0:jpl)
	761	!
[9421]	762	IF( hi_max(1) < rn_himin ) CALL ctl_stop('ice_itd_init: the upper bound of the 1st category must be bigger than rn_himin')
	763	!
[8586]	764	END SUBROUTINE ice_itd_init
	765
	766	#else
	767	!!----------------------------------------------------------------------
[9570]	768	!! Default option : Empty module NO SI3 sea-ice model
[8586]	769	!!----------------------------------------------------------------------
	770	#endif
	771
	772	!!======================================================================
	773	END MODULE iceitd

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

Download in other formats: