Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

limitd_th.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 33.5 KB

Rev	Line
[825]	1	MODULE limitd_th
	2	!!======================================================================
	3	!! * MODULE limitd_th *
[3625]	4	!! LIM3 ice model : ice thickness distribution: Thermodynamics
[825]	5	!!======================================================================
[2715]	6	!! History : - ! (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
	7	!! 3.0 ! 2005-12 (M. Vancoppenolle) adaptation to LIM-3
[4869]	8	!! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age
[2715]	9	!! - ! 2007-04 (M. Vancoppenolle) Mass conservation checked
	10	!!----------------------------------------------------------------------
[2528]	11	#if defined key_lim3
[825]	12	!!----------------------------------------------------------------------
[2528]	13	!! 'key_lim3' : LIM3 sea-ice model
	14	!!----------------------------------------------------------------------
[2715]	15	!! lim_itd_th_rem :
	16	!! lim_itd_th_reb :
	17	!! lim_itd_fitline :
	18	!! lim_itd_shiftice :
[2528]	19	!!----------------------------------------------------------------------
[4161]	20	USE par_oce ! ocean parameters
	21	USE dom_oce ! ocean domain
	22	USE phycst ! physical constants (ocean directory)
	23	USE thd_ice ! LIM-3 thermodynamic variables
	24	USE ice ! LIM-3 variables
	25	USE limvar ! LIM-3 variables
	26	USE prtctl ! Print control
	27	USE in_out_manager ! I/O manager
	28	USE lib_mpp ! MPP library
	29	USE lib_fortran ! to use key_nosignedzero
[5123]	30	USE limcons ! conservation tests
[921]	31
[825]	32	IMPLICIT NONE
	33	PRIVATE
	34
[2715]	35	PUBLIC lim_itd_th_rem
	36	PUBLIC lim_itd_th_reb
[825]	37
	38	!!----------------------------------------------------------------------
[4161]	39	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
[1156]	40	!! $Id$
[2715]	41	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	42	!!----------------------------------------------------------------------
	43	CONTAINS
	44
[4869]	45	SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
[921]	46	!!------------------------------------------------------------------
	47	!! * ROUTINE lim_itd_th_rem *
	48	!!
[2715]	49	!! ** Purpose : computes the redistribution of ice thickness
	50	!! after thermodynamic growth of ice thickness
	51	!!
[921]	52	!! ** Method : Linear remapping
	53	!!
[2715]	54	!! References : W.H. Lipscomb, JGR 2001
[921]	55	!!------------------------------------------------------------------
[2715]	56	INTEGER , INTENT (in) :: klbnd ! Start thickness category index point
	57	INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied
	58	INTEGER , INTENT (in) :: kt ! Ocean time step
	59	!
	60	INTEGER :: ji, jj, jl ! dummy loop index
[4161]	61	INTEGER :: ii, ij ! 2D corresponding indices to ji
	62	INTEGER :: nd ! local integer
[2715]	63	REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars
[4161]	64	REAL(wp) :: zx2, zwk2, zda0, zetamax ! - -
[5123]	65	REAL(wp) :: zx3
[2715]	66	CHARACTER (len = 15) :: fieldid
[825]	67
[7910]	68	INTEGER , DIMENSION(jpi,jpj,jpl-1) :: zdonor ! donor category index
[825]	69
[7910]	70	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zdhice ! ice thickness increment
	71	REAL(wp), DIMENSION(jpi,jpj,jpl) :: g0 ! coefficients for fitting the line of the ITD
	72	REAL(wp), DIMENSION(jpi,jpj,jpl) :: g1 ! coefficients for fitting the line of the ITD
	73	REAL(wp), DIMENSION(jpi,jpj,jpl) :: hL ! left boundary for the ITD for each thickness
	74	REAL(wp), DIMENSION(jpi,jpj,jpl) :: hR ! left boundary for the ITD for each thickness
	75	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zht_i_b ! old ice thickness
	76	REAL(wp), DIMENSION(jpi,jpj,jpl) :: dummy_es
	77	REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: zdaice, zdvice ! local increment of ice area and volume
	78	REAL(wp), DIMENSION((jpi+1)*(jpj+1)) :: zvetamin, zvetamax ! maximum values for etas
	79	INTEGER , DIMENSION((jpi+1)*(jpj+1)) :: nind_i, nind_j ! compressed indices for i/j directions
[3294]	80	INTEGER :: nbrem ! number of cells with ice to transfer
	81	REAL(wp) :: zslope ! used to compute local thermodynamic "speeds"
[7910]	82	REAL(wp), DIMENSION(jpi,jpj) :: zhb0, zhb1 ! category boundaries for thinnes categories
	83	REAL(wp), DIMENSION(jpi,jpj) :: vt_i_init, vt_i_final ! ice volume summed over categories
	84	REAL(wp), DIMENSION(jpi,jpj) :: vt_s_init, vt_s_final ! snow volume summed over categories
	85	REAL(wp), DIMENSION(jpi,jpj) :: et_i_init, et_i_final ! ice energy summed over categories
	86	REAL(wp), DIMENSION(jpi,jpj) :: et_s_init, et_s_final ! snow energy summed over categories
	87	INTEGER , DIMENSION(jpi,jpj) :: zremap_flag ! compute remapping or not ????
	88	REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: zhbnew ! new boundaries of ice categories
[3294]	89	!!------------------------------------------------------------------
[825]	90
	91
[921]	92	!!----------------------------------------------------------------------------------------------
	93	!! 0) Conservation checkand changes in each ice category
	94	!!----------------------------------------------------------------------------------------------
[2715]	95	IF( con_i ) THEN
[825]	96	CALL lim_column_sum (jpl, v_i, vt_i_init)
	97	CALL lim_column_sum (jpl, v_s, vt_s_init)
	98	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
[7753]	99	dummy_es(:,:,:) = e_s(:,:,1,:)
[825]	100	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
	101	ENDIF
	102
[921]	103	!!----------------------------------------------------------------------------------------------
	104	!! 1) Compute thickness and changes in each ice category
	105	!!----------------------------------------------------------------------------------------------
[2715]	106	IF( kt == nit000 .AND. lwp) THEN
[921]	107	WRITE(numout,*)
	108	WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution'
	109	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	110	WRITE(numout,*) ' klbnd : ', klbnd
	111	WRITE(numout,*) ' kubnd : ', kubnd
	112	ENDIF
[825]	113
[7753]	114	zdhice(:,:,:) = 0._wp
[921]	115	DO jl = klbnd, kubnd
	116	DO jj = 1, jpj
	117	DO ji = 1, jpi
[5134]	118	rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes
[4990]	119	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
[5407]	120	rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) )
[4990]	121	zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
[5407]	122	IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement?
[921]	123	END DO
	124	END DO
	125	END DO
	126
	127	!-----------------------------------------------------------------------------------------------
	128	! 2) Compute fractional ice area in each grid cell
	129	!-----------------------------------------------------------------------------------------------
[7753]	130	at_i(:,:) = 0._wp
[825]	131	DO jl = klbnd, kubnd
[7753]	132	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
[825]	133	END DO
	134
[921]	135	!-----------------------------------------------------------------------------------------------
	136	! 3) Identify grid cells with ice
	137	!-----------------------------------------------------------------------------------------------
[825]	138	nbrem = 0
	139	DO jj = 1, jpj
	140	DO ji = 1, jpi
[5123]	141	IF ( at_i(ji,jj) > epsi10 ) THEN
[825]	142	nbrem = nbrem + 1
	143	nind_i(nbrem) = ji
	144	nind_j(nbrem) = jj
[3294]	145	zremap_flag(ji,jj) = 1
[825]	146	ELSE
[3294]	147	zremap_flag(ji,jj) = 0
[825]	148	ENDIF
[5123]	149	END DO
	150	END DO
[825]	151
[921]	152	!-----------------------------------------------------------------------------------------------
	153	! 4) Compute new category boundaries
	154	!-----------------------------------------------------------------------------------------------
[825]	155	!- 4.1 Compute category boundaries
[7753]	156	zhbnew(:,:,:) = 0._wp
[825]	157
	158	DO jl = klbnd, kubnd - 1
	159	DO ji = 1, nbrem
[4161]	160	ii = nind_i(ji)
	161	ij = nind_j(ji)
[825]	162	!
[4688]	163	zhbnew(ii,ij,jl) = hi_max(jl)
[5407]	164	IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
[825]	165	!interpolate between adjacent category growth rates
[4872]	166	zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
	167	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
[5407]	168	ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN
[4161]	169	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
[5407]	170	ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN
[4161]	171	zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
[825]	172	ENDIF
[2715]	173	END DO
[825]	174
[921]	175	!- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
[825]	176	DO ji = 1, nbrem
[4161]	177	ii = nind_i(ji)
	178	ij = nind_j(ji)
[5407]	179
	180	! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible
	181	! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
	182	IF ( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN
[4161]	183	zremap_flag(ii,ij) = 0
[5407]	184	ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN
[4161]	185	zremap_flag(ii,ij) = 0
[825]	186	ENDIF
	187
[921]	188	!- 4.3 Check that each zhbnew does not exceed maximal values hi_max
[5407]	189	IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
[4688]	190	IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0
[5407]	191	! clem bug: why is not the following instead?
	192	!!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0
	193	!!IF( zhbnew(ii,ij,jl) > hi_max(jl ) ) zremap_flag(ii,ij) = 0
	194
[4688]	195	END DO
	196
[5123]	197	END DO
[825]	198
[921]	199	!-----------------------------------------------------------------------------------------------
	200	! 5) Identify cells where ITD is to be remapped
	201	!-----------------------------------------------------------------------------------------------
	202	nbrem = 0
	203	DO jj = 1, jpj
	204	DO ji = 1, jpi
[4688]	205	IF( zremap_flag(ji,jj) == 1 ) THEN
[921]	206	nbrem = nbrem + 1
	207	nind_i(nbrem) = ji
	208	nind_j(nbrem) = jj
	209	ENDIF
[4688]	210	END DO
	211	END DO
[825]	212
[921]	213	!-----------------------------------------------------------------------------------------------
	214	! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
	215	!-----------------------------------------------------------------------------------------------
	216	DO jj = 1, jpj
	217	DO ji = 1, jpi
[5407]	218	zhb0(ji,jj) = hi_max(0)
	219	zhb1(ji,jj) = hi_max(1)
[825]	220
[2715]	221	IF( a_i(ji,jj,kubnd) > epsi10 ) THEN
[5134]	222	zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) )
[921]	223	ELSE
[5407]	224	!clem bug zhbnew(ji,jj,kubnd) = hi_max(kubnd)
	225	zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway
[921]	226	ENDIF
[825]	227
[5407]	228	! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible
	229	! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
	230	IF ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) ) THEN
	231	zremap_flag(ji,jj) = 0
	232	ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) ) THEN
	233	zremap_flag(ji,jj) = 0
	234	ENDIF
	235
[5134]	236	END DO
	237	END DO
[825]	238
[921]	239	!-----------------------------------------------------------------------------------------------
	240	! 7) Compute g(h)
	241	!-----------------------------------------------------------------------------------------------
	242	!- 7.1 g(h) for category 1 at start of time step
[5123]	243	CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &
[2715]	244	& hR(:,:,klbnd), zremap_flag )
[825]	245
[921]	246	!- 7.2 Area lost due to melting of thin ice (first category, klbnd)
	247	DO ji = 1, nbrem
[4161]	248	ii = nind_i(ji)
	249	ij = nind_j(ji)
[825]	250
[5123]	251	IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
[5407]	252
[4161]	253	zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
[5407]	254
	255	IF( zdh0 < 0.0 ) THEN !remove area from category 1
[5123]	256	zdh0 = MIN( -zdh0, hi_max(klbnd) )
[921]	257	!Integrate g(1) from 0 to dh0 to estimate area melted
[5123]	258	zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
[5407]	259
[5123]	260	IF( zetamax > 0.0 ) THEN
[5407]	261	zx1 = zetamax
	262	zx2 = 0.5 * zetamax * zetamax
	263	zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 ! ice area removed
	264	zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i
	265	zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting
	266	! of thin ice (zdamax > 0)
[921]	267	! Remove area, conserving volume
[5123]	268	ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
[4161]	269	a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0
[5123]	270	v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
	271	ENDIF
[825]	272
[5407]	273	ELSE ! if ice accretion zdh0 > 0
	274	! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1)
[5123]	275	zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )
[5407]	276	ENDIF
[825]	277
[5407]	278	ENDIF
[825]	279
[5134]	280	END DO
[825]	281
[921]	282	!- 7.3 g(h) for each thickness category
	283	DO jl = klbnd, kubnd
[5123]	284	CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &
	285	& g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
[921]	286	END DO
[825]	287
[921]	288	!-----------------------------------------------------------------------------------------------
	289	! 8) Compute area and volume to be shifted across each boundary
	290	!-----------------------------------------------------------------------------------------------
[825]	291
[921]	292	DO jl = klbnd, kubnd - 1
	293	DO jj = 1, jpj
	294	DO ji = 1, jpi
	295	zdonor(ji,jj,jl) = 0
	296	zdaice(ji,jj,jl) = 0.0
	297	zdvice(ji,jj,jl) = 0.0
	298	END DO
	299	END DO
[825]	300
[921]	301	DO ji = 1, nbrem
[4161]	302	ii = nind_i(ji)
	303	ij = nind_j(ji)
[825]	304
[5123]	305	IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
[921]	306	! left and right integration limits in eta space
[5123]	307	zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
[5407]	308	zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
[4161]	309	zdonor(ii,ij,jl) = jl
[825]	310
[5407]	311	ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
[921]	312	! left and right integration limits in eta space
	313	zvetamin(ji) = 0.0
[5123]	314	zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
[4161]	315	zdonor(ii,ij,jl) = jl + 1
[825]	316
[5407]	317	ENDIF
[825]	318
[5123]	319	zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
[921]	320	zetamin = zvetamin(ji)
[825]	321
[921]	322	zx1 = zetamax - zetamin
[5123]	323	zwk1 = zetamin * zetamin
	324	zwk2 = zetamax * zetamax
	325	zx2 = 0.5 * ( zwk2 - zwk1 )
[921]	326	zwk1 = zwk1 * zetamin
	327	zwk2 = zwk2 * zetamax
[5123]	328	zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 )
[4161]	329	nd = zdonor(ii,ij,jl)
	330	zdaice(ii,ij,jl) = g1(ii,ij,nd)zx2 + g0(ii,ij,nd)zx1
[4688]	331	zdvice(ii,ij,jl) = g1(ii,ij,nd)zx3 + g0(ii,ij,nd)zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
[921]	332
[5123]	333	END DO
[5407]	334	END DO
[921]	335
	336	!!----------------------------------------------------------------------------------------------
	337	!! 9) Shift ice between categories
	338	!!----------------------------------------------------------------------------------------------
	339	CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
	340
	341	!!----------------------------------------------------------------------------------------------
	342	!! 10) Make sure ht_i >= minimum ice thickness hi_min
	343	!!----------------------------------------------------------------------------------------------
	344
	345	DO ji = 1, nbrem
[4161]	346	ii = nind_i(ji)
	347	ij = nind_j(ji)
[5123]	348	IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
[5134]	349	a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin
[5123]	350	ht_i(ii,ij,1) = rn_himin
[921]	351	ENDIF
[5123]	352	END DO
[921]	353
	354	!!----------------------------------------------------------------------------------------------
	355	!! 11) Conservation check
	356	!!----------------------------------------------------------------------------------------------
[825]	357	IF ( con_i ) THEN
	358	CALL lim_column_sum (jpl, v_i, vt_i_final)
	359	fieldid = ' v_i : limitd_th '
	360	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
	361
	362	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_final)
	363	fieldid = ' e_i : limitd_th '
	364	CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid)
	365
	366	CALL lim_column_sum (jpl, v_s, vt_s_final)
	367	fieldid = ' v_s : limitd_th '
	368	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
	369
[7753]	370	dummy_es(:,:,:) = e_s(:,:,1,:)
[825]	371	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
	372	fieldid = ' e_s : limitd_th '
	373	CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
	374	ENDIF
	375
[3294]	376
[921]	377	END SUBROUTINE lim_itd_th_rem
[825]	378
	379
[5123]	380	SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
[921]	381	!!------------------------------------------------------------------
	382	!! * ROUTINE lim_itd_fitline *
	383	!!
[2715]	384	!! ** Purpose : fit g(h) with a line using area, volume constraints
[921]	385	!!
[2715]	386	!! ** Method : Fit g(h) with a line, satisfying area and volume constraints.
	387	!! To reduce roundoff errors caused by large values of g0 and g1,
	388	!! we actually compute g(eta), where eta = h - hL, and hL is the
	389	!! left boundary.
[921]	390	!!------------------------------------------------------------------
[2715]	391	INTEGER , INTENT(in ) :: num_cat ! category index
	392	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: HbL, HbR ! left and right category boundaries
	393	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hice ! ice thickness
	394	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: g0, g1 ! coefficients in linear equation for g(eta)
	395	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hL ! min value of range over which g(h) > 0
	396	REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hR ! max value of range over which g(h) > 0
[3294]	397	INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag !
[2715]	398	!
[5407]	399	INTEGER :: ji,jj ! horizontal indices
[2715]	400	REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL)
	401	REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL)
	402	REAL(wp) :: zdhr ! 1 / (hR - hL)
	403	REAL(wp) :: zwk1, zwk2 ! temporary variables
	404	!!------------------------------------------------------------------
	405	!
[825]	406	DO jj = 1, jpj
	407	DO ji = 1, jpi
[2715]	408	!
[4333]	409	IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 &
[5407]	410	& .AND. hice(ji,jj) > 0._wp ) THEN
[825]	411
[921]	412	! Initialize hL and hR
[825]	413	hL(ji,jj) = HbL(ji,jj)
	414	hR(ji,jj) = HbR(ji,jj)
	415
[921]	416	! Change hL or hR if hice falls outside central third of range
[5123]	417	zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
	418	zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
[825]	419
[2715]	420	IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
	421	ELSEIF( hice(ji,jj) > zh23 ) THEN ; hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj)
[825]	422	ENDIF
	423
[921]	424	! Compute coefficients of g(eta) = g0 + g1*eta
[2715]	425	zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj))
	426	zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
	427	zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
[5123]	428	g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
	429	g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
[2715]	430	!
[5407]	431	ELSE ! remap_flag = .false. or a_i < epsi10
[2715]	432	hL(ji,jj) = 0._wp
	433	hR(ji,jj) = 0._wp
	434	g0(ji,jj) = 0._wp
	435	g1(ji,jj) = 0._wp
[5407]	436	ENDIF
[2715]	437	!
	438	END DO
	439	END DO
	440	!
	441	END SUBROUTINE lim_itd_fitline
[825]	442
	443
[2715]	444	SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
[921]	445	!!------------------------------------------------------------------
	446	!! * ROUTINE lim_itd_shiftice *
[2715]	447	!!
	448	!! ** Purpose : shift ice across category boundaries, conserving everything
[921]	449	!! ( area, volume, energy, age*vol, and mass of salt )
	450	!!
	451	!! ** Method :
	452	!!------------------------------------------------------------------
[3294]	453	INTEGER , INTENT(in ) :: klbnd ! Start thickness category index point
	454	INTEGER , INTENT(in ) :: kubnd ! End point on which the the computation is applied
[2715]	455	INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in ) :: zdonor ! donor category index
	456	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdaice ! ice area transferred across boundary
	457	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdvice ! ice volume transferred across boundary
[825]	458
[2715]	459	INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices
[5407]	460	INTEGER :: ii, ij ! indices when changing from 2D-1D is done
[825]	461
[7910]	462	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zaTsfn
	463	REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here
[825]	464
[2715]	465	REAL(wp) :: zdvsnow, zdesnow ! snow volume and energy transferred
	466	REAL(wp) :: zdeice ! ice energy transferred
	467	REAL(wp) :: zdsm_vice ! ice salinity times volume transferred
	468	REAL(wp) :: zdo_aice ! ice age times volume transferred
	469	REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred
[825]	470
[7910]	471	INTEGER, DIMENSION((jpi+1)*(jpj+1)) :: nind_i, nind_j ! compressed indices for i/j directions
[825]	472
[5407]	473	INTEGER :: nbrem ! number of cells with ice to transfer
[2715]	474	!!------------------------------------------------------------------
[825]	475
[3294]	476
[921]	477	!----------------------------------------------------------------------------------------------
	478	! 1) Define a variable equal to a_i*T_su
	479	!----------------------------------------------------------------------------------------------
[825]	480
	481	DO jl = klbnd, kubnd
[7753]	482	zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
[2715]	483	END DO
[825]	484
[921]	485	!-------------------------------------------------------------------------------
[5407]	486	! 2) Transfer volume and energy between categories
[921]	487	!-------------------------------------------------------------------------------
[825]	488
	489	DO jl = klbnd, kubnd - 1
	490	nbrem = 0
	491	DO jj = 1, jpj
	492	DO ji = 1, jpi
[5123]	493	IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
[825]	494	nbrem = nbrem + 1
	495	nind_i(nbrem) = ji
	496	nind_j(nbrem) = jj
[5123]	497	ENDIF
[825]	498	END DO
	499	END DO
	500
	501	DO ji = 1, nbrem
[4161]	502	ii = nind_i(ji)
	503	ij = nind_j(ji)
[825]	504
[4161]	505	jl1 = zdonor(ii,ij,jl)
[5407]	506	rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) )
	507	zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch
[2715]	508	IF( jl1 == jl) THEN ; jl2 = jl1+1
[5123]	509	ELSE ; jl2 = jl
[825]	510	ENDIF
	511
	512	!--------------
	513	! Ice areas
	514	!--------------
[4161]	515	a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl)
	516	a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl)
[825]	517
	518	!--------------
	519	! Ice volumes
	520	!--------------
[4161]	521	v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)
	522	v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl)
[825]	523
	524	!--------------
	525	! Snow volumes
	526	!--------------
[4688]	527	zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij)
[4161]	528	v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow
	529	v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow
[825]	530
	531	!--------------------
	532	! Snow heat content
	533	!--------------------
[4688]	534	zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij)
[4161]	535	e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow
	536	e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow
[825]	537
	538	!--------------
	539	! Ice age
	540	!--------------
[4688]	541	zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl)
[4161]	542	oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice
	543	oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice
[825]	544
	545	!--------------
	546	! Ice salinity
	547	!--------------
[4688]	548	zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij)
[4161]	549	smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice
	550	smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice
[825]	551
	552	!---------------------
	553	! Surface temperature
	554	!---------------------
[4688]	555	zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl)
[4161]	556	zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf
	557	zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf
[825]	558
[5123]	559	END DO
[825]	560
	561	!------------------
	562	! Ice heat content
	563	!------------------
	564
	565	DO jk = 1, nlay_i
	566	DO ji = 1, nbrem
[4161]	567	ii = nind_i(ji)
	568	ij = nind_j(ji)
[825]	569
[4161]	570	jl1 = zdonor(ii,ij,jl)
[5123]	571	IF (jl1 == jl) THEN
[825]	572	jl2 = jl+1
	573	ELSE ! n1 = n+1
	574	jl2 = jl
	575	ENDIF
	576
[4161]	577	zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij)
	578	e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice
	579	e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice
[5123]	580	END DO
	581	END DO
[825]	582
	583	END DO ! boundaries, 1 to ncat-1
	584
	585	!-----------------------------------------------------------------
	586	! Update ice thickness and temperature
	587	!-----------------------------------------------------------------
	588
	589	DO jl = klbnd, kubnd
	590	DO jj = 1, jpj
[921]	591	DO ji = 1, jpi
[2715]	592	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
	593	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
[921]	594	t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)
	595	ELSE
[2715]	596	ht_i(ji,jj,jl) = 0._wp
[5123]	597	t_su(ji,jj,jl) = rt0
[921]	598	ENDIF
[5123]	599	END DO
	600	END DO
	601	END DO
[2715]	602	!
[3294]	603	!
[921]	604	END SUBROUTINE lim_itd_shiftice
[2715]	605
[825]	606
[4869]	607	SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
[921]	608	!!------------------------------------------------------------------
	609	!! * ROUTINE lim_itd_th_reb *
[2715]	610	!!
[921]	611	!! ** Purpose : rebin - rebins thicknesses into defined categories
	612	!!
	613	!! ** Method :
	614	!!------------------------------------------------------------------
[2715]	615	INTEGER , INTENT (in) :: klbnd ! Start thickness category index point
	616	INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied
	617	!
	618	INTEGER :: ji,jj, jl ! dummy loop indices
	619	INTEGER :: zshiftflag ! = .true. if ice must be shifted
[921]	620	CHARACTER (len = 15) :: fieldid
[825]	621
[7910]	622	INTEGER , DIMENSION(jpi,jpj,jpl) :: zdonor ! donor category index
	623	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zdaice, zdvice ! ice area and volume transferred
[825]	624
[7910]	625	REAL(wp), DIMENSION(jpi,jpj) :: vt_i_init, vt_i_final ! ice volume summed over categories
	626	REAL(wp), DIMENSION(jpi,jpj) :: vt_s_init, vt_s_final ! snow volume summed over categories
[2715]	627	!!------------------------------------------------------------------
[3294]	628
[2715]	629	!
	630	IF( con_i ) THEN ! conservation check
[834]	631	CALL lim_column_sum (jpl, v_i, vt_i_init)
	632	CALL lim_column_sum (jpl, v_s, vt_s_init)
	633	ENDIF
[825]	634
[921]	635	!
	636	!------------------------------------------------------------------------------
	637	! 1) Compute ice thickness.
	638	!------------------------------------------------------------------------------
[825]	639	DO jl = klbnd, kubnd
	640	DO jj = 1, jpj
[921]	641	DO ji = 1, jpi
[5134]	642	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
	643	ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch
[2715]	644	END DO
	645	END DO
	646	END DO
[825]	647
[921]	648	!------------------------------------------------------------------------------
[5134]	649	! 2) If a category thickness is not in bounds, shift the
[921]	650	! entire area, volume, and energy to the neighboring category
	651	!------------------------------------------------------------------------------
[825]	652	!-------------------------
	653	! Initialize shift arrays
	654	!-------------------------
	655	DO jl = klbnd, kubnd
[7753]	656	zdonor(:,:,jl) = 0
	657	zdaice(:,:,jl) = 0._wp
	658	zdvice(:,:,jl) = 0._wp
[825]	659	END DO
	660
	661	!-------------------------
	662	! Move thin categories up
	663	!-------------------------
	664
	665	DO jl = klbnd, kubnd - 1 ! loop over category boundaries
	666
[921]	667	!---------------------------------------
	668	! identify thicknesses that are too big
	669	!---------------------------------------
[869]	670	zshiftflag = 0
[825]	671
	672	DO jj = 1, jpj
	673	DO ji = 1, jpi
[2715]	674	IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN
[869]	675	zshiftflag = 1
[825]	676	zdonor(ji,jj,jl) = jl
[4161]	677	! begin TECLIM change
[4293]	678	!zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp
	679	!zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)(hi_max(jl)+hi_max(jl-1)) 0.5_wp
[4161]	680	! end TECLIM change
[4293]	681	! clem: how much of a_i you send in cat sup is somewhat arbitrary
[5134]	682	zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl)
	683	zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 )
[825]	684	ENDIF
[5123]	685	END DO
	686	END DO
[2715]	687	IF(lk_mpp) CALL mpp_max( zshiftflag )
[825]	688
[2715]	689	IF( zshiftflag == 1 ) THEN ! Shift ice between categories
	690	CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
[921]	691	! Reset shift parameters
[7753]	692	zdonor(:,:,jl) = 0
	693	zdaice(:,:,jl) = 0._wp
	694	zdvice(:,:,jl) = 0._wp
[2715]	695	ENDIF
	696	!
[5134]	697	END DO
[825]	698
	699	!----------------------------
	700	! Move thick categories down
	701	!----------------------------
	702
	703	DO jl = kubnd - 1, 1, -1 ! loop over category boundaries
	704
[921]	705	!-----------------------------------------
	706	! Identify thicknesses that are too small
	707	!-----------------------------------------
[869]	708	zshiftflag = 0
[825]	709
[4688]	710	DO jj = 1, jpj
	711	DO ji = 1, jpi
	712	IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN
[7753]	713	!
[4688]	714	zshiftflag = 1
	715	zdonor(ji,jj,jl) = jl + 1
	716	zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)
	717	zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
	718	ENDIF
[5123]	719	END DO
	720	END DO
[4688]	721
	722	IF(lk_mpp) CALL mpp_max( zshiftflag )
	723
	724	IF( zshiftflag == 1 ) THEN ! Shift ice between categories
	725	CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )
	726	! Reset shift parameters
[7753]	727	zdonor(:,:,jl) = 0
	728	zdaice(:,:,jl) = 0._wp
	729	zdvice(:,:,jl) = 0._wp
[4688]	730	ENDIF
	731
[5123]	732	END DO
[825]	733
[921]	734	!------------------------------------------------------------------------------
[5134]	735	! 3) Conservation check
[921]	736	!------------------------------------------------------------------------------
[825]	737
[2715]	738	IF( con_i ) THEN
[921]	739	CALL lim_column_sum (jpl, v_i, vt_i_final)
	740	fieldid = ' v_i : limitd_reb '
	741	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
[825]	742
[921]	743	CALL lim_column_sum (jpl, v_s, vt_s_final)
	744	fieldid = ' v_s : limitd_reb '
	745	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
	746	ENDIF
[2715]	747	!
[3294]	748
[921]	749	END SUBROUTINE lim_itd_th_reb
[825]	750
	751	#else
[2715]	752	!!----------------------------------------------------------------------
	753	!! Default option Dummy module NO LIM sea-ice model
	754	!!----------------------------------------------------------------------
[825]	755	CONTAINS
	756	SUBROUTINE lim_itd_th_rem
	757	END SUBROUTINE lim_itd_th_rem
	758	SUBROUTINE lim_itd_fitline
	759	END SUBROUTINE lim_itd_fitline
	760	SUBROUTINE lim_itd_shiftice
	761	END SUBROUTINE lim_itd_shiftice
	762	SUBROUTINE lim_itd_th_reb
	763	END SUBROUTINE lim_itd_th_reb
	764	#endif
[2715]	765	!!======================================================================
[921]	766	END MODULE limitd_th

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90 @ 7910

Download in other formats: