Context Navigation

limitd_th.F90 @ 7753

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

Download in other formats: