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.

limthd_lac.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 4333

Last change on this file since 4333 was 4333, checked in by clem, 10 years ago
remove remaining bugs in LIM3, so that it can run in both regional and global config
Property svn:keywords set to `Id`
File size: 33.2 KB

Rev	Line
[825]	1	MODULE limthd_lac
	2	!!======================================================================
	3	!! * MODULE limthd_lac *
	4	!! lateral thermodynamic growth of the ice
	5	!!======================================================================
[2715]	6	!! History : LIM ! 2005-12 (M. Vancoppenolle) Original code
	7	!! - ! 2006-01 (M. Vancoppenolle) add ITD
	8	!! 3.0 ! 2007-07 (M. Vancoppenolle) Mass and energy conservation tested
	9	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	10	!!----------------------------------------------------------------------
[888]	11	#if defined key_lim3
[825]	12	!!----------------------------------------------------------------------
[2528]	13	!! 'key_lim3' LIM3 sea-ice model
	14	!!----------------------------------------------------------------------
[3625]	15	!! lim_lat_acr : lateral accretion of ice
[2528]	16	!!----------------------------------------------------------------------
[3625]	17	USE par_oce ! ocean parameters
	18	USE dom_oce ! domain variables
	19	USE phycst ! physical constants
	20	USE sbc_oce ! Surface boundary condition: ocean fields
	21	USE sbc_ice ! Surface boundary condition: ice fields
	22	USE thd_ice ! LIM thermodynamics
	23	USE dom_ice ! LIM domain
	24	USE par_ice ! LIM parameters
	25	USE ice ! LIM variables
	26	USE limtab ! LIM 2D <==> 1D
	27	USE limcons ! LIM conservation
	28	USE in_out_manager ! I/O manager
	29	USE lib_mpp ! MPP library
	30	USE wrk_nemo ! work arrays
	31	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[921]	32
[825]	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC lim_thd_lac ! called by lim_thd
	37
[4333]	38	REAL(wp) :: epsi10 = 1.e-10_wp !
[2715]	39	REAL(wp) :: zzero = 0._wp !
	40	REAL(wp) :: zone = 1._wp !
[825]	41
	42	!!----------------------------------------------------------------------
[4161]	43	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	44	!! $Id$
[2715]	45	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	46	!!----------------------------------------------------------------------
	47	CONTAINS
[921]	48
[825]	49	SUBROUTINE lim_thd_lac
	50	!!-------------------------------------------------------------------
	51	!! * ROUTINE lim_thd_lac *
	52	!!
	53	!! ** Purpose : Computation of the evolution of the ice thickness and
	54	!! concentration as a function of the heat balance in the leads.
	55	!! It is only used for lateral accretion
	56	!!
	57	!! ** Method : Ice is formed in the open water when ocean lose heat
	58	!! (heat budget of open water Bl is negative) .
	59	!! Computation of the increase of 1-A (ice concentration) fol-
	60	!! lowing the law :
	61	!! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
	62	!! where - h0 is the thickness of ice created in the lead
	63	!! - a is a minimum fraction for leads
	64	!! - F is a monotonic non-increasing function defined as:
	65	!! F(X)=( 1 - Xexld )(1.0/exld)
	66	!! - exld is the exponent closure rate (=2 default val.)
	67	!!
	68	!! ** Action : - Adjustment of snow and ice thicknesses and heat
	69	!! content in brine pockets
	70	!! - Updating ice internal temperature
	71	!! - Computation of variation of ice volume and mass
	72	!! - Computation of frldb after lateral accretion and
	73	!! update ht_s_b, ht_i_b and tbif_1d(:,:)
	74	!!------------------------------------------------------------------------
[4161]	75	INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices
	76	INTEGER :: layer, nbpac ! local integers
	77	INTEGER :: ii, ij, iter ! - -
	78	REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars
[2715]	79	REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - -
	80	REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - -
	81	LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness
	82	CHARACTER (len = 15) :: fieldid
	83	!
[3294]	84	INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows
	85	REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not
[825]	86
[3294]	87	REAL(wp), POINTER, DIMENSION(:) :: zv_newice ! volume of accreted ice
	88	REAL(wp), POINTER, DIMENSION(:) :: za_newice ! fractional area of accreted ice
	89	REAL(wp), POINTER, DIMENSION(:) :: zh_newice ! thickness of accreted ice
	90	REAL(wp), POINTER, DIMENSION(:) :: ze_newice ! heat content of accreted ice
	91	REAL(wp), POINTER, DIMENSION(:) :: zs_newice ! salinity of accreted ice
	92	REAL(wp), POINTER, DIMENSION(:) :: zo_newice ! age of accreted ice
	93	REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget
	94	REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget
	95	REAL(wp), POINTER, DIMENSION(:) :: zat_i_ac ! total ice fraction
	96	REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1)
	97	REAL(wp), POINTER, DIMENSION(:) :: zdh_frazb ! accretion of frazil ice at the ice bottom
	98	REAL(wp), POINTER, DIMENSION(:) :: zvrel_ac ! relative ice / frazil velocity (1D vector)
[825]	99
[3294]	100	REAL(wp), POINTER, DIMENSION(:,:) :: zhice_old ! previous ice thickness
	101	REAL(wp), POINTER, DIMENSION(:,:) :: zdummy ! dummy thickness of new ice
	102	REAL(wp), POINTER, DIMENSION(:,:) :: zdhicbot ! thickness of new ice which is accreted vertically
	103	REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl
	104	REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl
	105	REAL(wp), POINTER, DIMENSION(:,:) :: za_i_ac ! 1-D version of a_i
	106	REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_ac ! 1-D version of v_i
	107	REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_ac ! 1-D version of oa_i
	108	REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_ac ! 1-D version of smv_i
[825]	109
[3294]	110	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_ac !: 1-D version of e_i
[825]	111
[3294]	112	REAL(wp), POINTER, DIMENSION(:) :: zqbgow ! heat budget of the open water (negative)
	113	REAL(wp), POINTER, DIMENSION(:) :: zdhex ! excessively thick accreted sea ice (hlead-hice)
[825]	114
[3294]	115	REAL(wp), POINTER, DIMENSION(:,:,:) :: zqm0 ! old layer-system heat content
	116	REAL(wp), POINTER, DIMENSION(:,:,:) :: zthick0 ! old ice thickness
[825]	117
[3294]	118	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
	119	REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories
	120	REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories
	121	REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init ! snow energy summed over categories
	122	REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity
	123	!!-----------------------------------------------------------------------!
[825]	124
[3294]	125	CALL wrk_alloc( jpij, zcatac ) ! integer
	126	CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
	127	CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
	128	CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
	129	CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
	130	CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
	131	CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
	132
[2715]	133	et_i_init(:,:) = 0._wp
	134	et_s_init(:,:) = 0._wp
	135	vt_i_init(:,:) = 0._wp
	136	vt_s_init(:,:) = 0._wp
[825]	137
[921]	138	!------------------------------------------------------------------------------!
	139	! 1) Conservation check and changes in each ice category
	140	!------------------------------------------------------------------------------!
[3625]	141	IF( con_i ) THEN
	142	CALL lim_column_sum ( jpl, v_i , vt_i_init)
	143	CALL lim_column_sum ( jpl, v_s , vt_s_init)
	144	CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
	145	CALL lim_column_sum ( jpl, e_s(:,:,1,:) , et_s_init)
[834]	146	ENDIF
[825]	147
[921]	148	!------------------------------------------------------------------------------\|
	149	! 2) Convert units for ice internal energy
	150	!------------------------------------------------------------------------------\|
[825]	151	DO jl = 1, jpl
[921]	152	DO jk = 1, nlay_i
	153	DO jj = 1, jpj
	154	DO ji = 1, jpi
	155	!Energy of melting q(S,T) [J.m-3]
[4161]	156	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i )
[4333]	157	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes
[3625]	158	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
[921]	159	END DO
[825]	160	END DO
[921]	161	END DO
[825]	162	END DO
	163
[921]	164	!------------------------------------------------------------------------------!
	165	! 3) Collection thickness of ice formed in leads and polynyas
	166	!------------------------------------------------------------------------------!
[865]	167	! hicol is the thickness of new ice formed in open water
	168	! hicol can be either prescribed (frazswi = 0)
	169	! or computed (frazswi = 1)
[825]	170	! Frazil ice forms in open water, is transported by wind
	171	! accumulates at the edge of the consolidated ice edge
	172	! where it forms aggregates of a specific thickness called
	173	! collection thickness.
	174
[865]	175	! Note : the following algorithm currently breaks vectorization
	176	!
	177
[3625]	178	zvrel(:,:) = 0._wp
[825]	179
	180	! Default new ice thickness
[3625]	181	hicol(:,:) = hiccrit(1)
[825]	182
[3625]	183	IF( fraz_swi == 1._wp ) THEN
[825]	184
[921]	185	!--------------------
	186	! Physical constants
	187	!--------------------
[3625]	188	hicol(:,:) = 0._wp
[825]	189
[921]	190	zhicrit = 0.04 ! frazil ice thickness
	191	ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
	192	zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
	193	zgamafr = 0.03
[825]	194
[921]	195	DO jj = 1, jpj
	196	DO ji = 1, jpi
[825]	197
[921]	198	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
	199	!-------------
	200	! Wind stress
	201	!-------------
	202	! C-grid wind stress components
[3625]	203	ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj ) &
	204	& + utau_ice(ji ,jj ) * tmu(ji ,jj ) ) * 0.5_wp
	205	ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &
	206	& + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp
[921]	207	! Square root of wind stress
	208	ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
[825]	209
[921]	210	!---------------------
	211	! Frazil ice velocity
	212	!---------------------
[2715]	213	zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
	214	zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
[825]	215
[921]	216	!-------------------
	217	! Pack ice velocity
	218	!-------------------
	219	! C-grid ice velocity
[3625]	220	zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) )
	221	zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) &
	222	& + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp
	223	zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) &
	224	& + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp
[825]	225
[921]	226	!-----------------------------------
	227	! Relative frazil/pack ice velocity
	228	!-----------------------------------
	229	! absolute relative velocity
[3625]	230	zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) &
	231	& + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
	232	zvrel(ji,jj) = SQRT( zvrel2 )
[825]	233
[921]	234	!---------------------
	235	! Iterative procedure
	236	!---------------------
	237	hicol(ji,jj) = zhicrit + 0.1
[3625]	238	hicol(ji,jj) = zhicrit + hicol(ji,jj) &
	239	& / ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) * ztwogp * zvrel2
[825]	240
[3625]	241	!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
	242	!!gm = zhicrit*2 + 0.2zhicrit +0.01
	243	!!gm therefore the 2 lines with hicol can be replaced by 1 line:
	244	!!gm hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
	245	!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
	246
[921]	247	iter = 1
	248	iterate_frazil = .true.
[825]	249
[921]	250	DO WHILE ( iter .LT. 100 .AND. iterate_frazil )
	251	zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)2 - zhicrit2 ) &
	252	- hicol(ji,jj) * zhicrit * ztwogp * zvrel2
	253	zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
	254	- zhicrit * ztwogp * zvrel2
	255	zhicol_new = hicol(ji,jj) - zf/zfp
	256	hicol(ji,jj) = zhicol_new
[825]	257
[921]	258	iter = iter + 1
[825]	259
[921]	260	END DO ! do while
[825]	261
[921]	262	ENDIF ! end of selection of pixels where ice forms
[825]	263
[921]	264	END DO ! loop on ji ends
	265	END DO ! loop on jj ends
[825]	266
	267	ENDIF ! End of computation of frazil ice collection thickness
	268
[921]	269	!------------------------------------------------------------------------------!
	270	! 4) Identify grid points where new ice forms
	271	!------------------------------------------------------------------------------!
[825]	272
	273	!-------------------------------------
	274	! Select points for new ice formation
	275	!-------------------------------------
	276	! This occurs if open water energy budget is negative
	277	nbpac = 0
	278	DO jj = 1, jpj
	279	DO ji = 1, jpi
[3625]	280	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN
[825]	281	nbpac = nbpac + 1
	282	npac( nbpac ) = (jj - 1) * jpi + ji
	283	ENDIF
	284	END DO
	285	END DO
	286
[4333]	287	! debug point to follow
	288	jiindex_1d = 0
	289	IF( ln_nicep ) THEN
	290	DO ji = mi0(jiindx), mi1(jiindx)
	291	DO jj = mj0(jjindx), mj1(jjindx)
	292	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN
	293	jiindex_1d = (jj - 1) * jpi + ji
	294	ENDIF
	295	END DO
	296	END DO
	297	ENDIF
	298
	299	IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
[825]	300
	301	!------------------------------
	302	! Move from 2-D to 1-D vectors
	303	!------------------------------
	304	! If ocean gains heat do nothing
	305	! 0therwise compute new ice formation
	306
	307	IF ( nbpac > 0 ) THEN
	308
[3625]	309	CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )
[921]	310	DO jl = 1, jpl
[3625]	311	CALL tab_2d_1d( nbpac, za_i_ac (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
	312	CALL tab_2d_1d( nbpac, zv_i_ac (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
	313	CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
	314	CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
[921]	315	DO jk = 1, nlay_i
[3625]	316	CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
[921]	317	END DO ! jk
	318	END DO ! jl
[825]	319
[3625]	320	CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )
	321	CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )
	322	CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )
	323	CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) )
	324	CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) )
	325	CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) )
	326	CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) )
[834]	327
[921]	328	!------------------------------------------------------------------------------!
	329	! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
	330	!------------------------------------------------------------------------------!
[825]	331
[921]	332	!----------------------
	333	! Thickness of new ice
	334	!----------------------
	335	DO ji = 1, nbpac
[3625]	336	zh_newice(ji) = hiccrit(1)
[921]	337	END DO
[3625]	338	IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:)
[825]	339
[921]	340	!----------------------
	341	! Salinity of new ice
	342	!----------------------
[825]	343
[3625]	344	SELECT CASE ( num_sal )
	345	CASE ( 1 ) ! Sice = constant
	346	zs_newice(:) = bulk_sal
	347	CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
[921]	348	DO ji = 1, nbpac
[4161]	349	ii = MOD( npac(ji) - 1 , jpi ) + 1
	350	ij = ( npac(ji) - 1 ) / jpi + 1
	351	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) )
[3625]	352	END DO
	353	CASE ( 3 ) ! Sice = F(z) [multiyear ice]
	354	zs_newice(:) = 2.3
	355	END SELECT
[825]	356
	357
[921]	358	!-------------------------
	359	! Heat content of new ice
	360	!-------------------------
	361	! We assume that new ice is formed at the seawater freezing point
	362	DO ji = 1, nbpac
[3625]	363	ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K)
	364	ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) &
	365	& + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) &
	366	& - rcp * ( ztmelts - rtt ) )
	367	ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) &
	368	& + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus
[921]	369	END DO ! ji
	370	!----------------
	371	! Age of new ice
	372	!----------------
	373	DO ji = 1, nbpac
[3625]	374	zo_newice(ji) = 0._wp
[921]	375	END DO ! ji
[825]	376
[921]	377	!--------------------------
	378	! Open water energy budget
	379	!--------------------------
	380	DO ji = 1, nbpac
[3625]	381	zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0
[921]	382	END DO ! ji
[825]	383
[921]	384	!-------------------
	385	! Volume of new ice
	386	!-------------------
	387	DO ji = 1, nbpac
[3625]	388	zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
[825]	389
[921]	390	! A fraction zfrazb of frazil ice is accreted at the ice bottom
[3625]	391	zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
	392	zdh_frazb(ji) = zfrazb * zv_newice(ji)
[921]	393	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
	394	END DO
[865]	395
[921]	396	!------------------------------------
	397	! Diags for energy conservation test
	398	!------------------------------------
	399	DO ji = 1, nbpac
[4161]	400	ii = MOD( npac(ji) - 1 , jpi ) + 1
	401	ij = ( npac(ji) - 1 ) / jpi + 1
[3625]	402	!
[4161]	403	zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
[3625]	404	!
[4161]	405	vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume
	406	et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy
[3625]	407
[921]	408	END DO
[825]	409
[921]	410	! keep new ice volume in memory
[3625]	411	CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
[825]	412
[921]	413	!-----------------
	414	! Area of new ice
	415	!-----------------
	416	DO ji = 1, nbpac
[4161]	417	ii = MOD( npac(ji) - 1 , jpi ) + 1
	418	ij = ( npac(ji) - 1 ) / jpi + 1
[3625]	419	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
[4161]	420	diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
[921]	421	END DO !ji
[825]	422
[921]	423	!------------------------------------------------------------------------------!
	424	! 6) Redistribute new ice area and volume into ice categories !
	425	!------------------------------------------------------------------------------!
[825]	426
[921]	427	!-----------------------------------------
	428	! Keep old ice areas and volume in memory
	429	!-----------------------------------------
	430	zv_old(:,:) = zv_i_ac(:,:)
	431	za_old(:,:) = za_i_ac(:,:)
[825]	432
[921]	433	!-------------------------------------------
	434	! Compute excessive new ice area and volume
	435	!-------------------------------------------
	436	! If lateral ice growth gives an ice concentration gt 1, then
[3625]	437	! we keep the excessive volume in memory and attribute it later to bottom accretion
[921]	438	DO ji = 1, nbpac
[4161]	439	IF ( za_newice(ji) > ( amax - zat_i_ac(ji) ) ) THEN
	440	zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) )
[3625]	441	zdv_res(ji) = zda_res (ji) * zh_newice(ji)
	442	za_newice(ji) = za_newice(ji) - zda_res (ji)
	443	zv_newice(ji) = zv_newice(ji) - zdv_res (ji)
[921]	444	ELSE
[3625]	445	zda_res(ji) = 0._wp
	446	zdv_res(ji) = 0._wp
[921]	447	ENDIF
	448	END DO ! ji
[825]	449
[921]	450	!------------------------------------------------
	451	! Laterally redistribute new ice volume and area
	452	!------------------------------------------------
[2715]	453	zat_i_ac(:) = 0._wp
[921]	454	DO jl = 1, jpl
	455	DO ji = 1, nbpac
[3625]	456	IF( hi_max (jl-1) < zh_newice(ji) .AND. &
	457	& zh_newice(ji) <= hi_max (jl) ) THEN
[2715]	458	za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
	459	zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
	460	zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl)
	461	zcatac (ji) = jl
[921]	462	ENDIF
[3625]	463	END DO
	464	END DO
[825]	465
[921]	466	!----------------------------------
	467	! Heat content - lateral accretion
	468	!----------------------------------
	469	DO ji = 1, nbpac
[2715]	470	jl = zcatac(ji) ! categroy in which new ice is put
[4161]	471	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise
[2715]	472	zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness
	473	zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness
[4161]	474	zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category
[921]	475	END DO
[825]	476
[921]	477	DO jk = 1, nlay_i
	478	DO ji = 1, nbpac
	479	jl = zcatac(ji)
[4161]	480	zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
	481	zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) &
	482	& - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
[2715]	483	ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) &
[4161]	484	+ ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) &
[2715]	485	+ za_newice(ji) * ze_newice(ji) * zalphai &
[4161]	486	+ za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
[2715]	487	END DO
	488	END DO
[825]	489
[921]	490	!-----------------------------------------------
	491	! Add excessive volume of new ice at the bottom
	492	!-----------------------------------------------
	493	! If the ice concentration exceeds 1, the remaining volume of new ice
	494	! is equally redistributed among all ice categories in which there is
	495	! ice
[825]	496
[921]	497	! Fraction of level ice
	498	jm = 1
[2715]	499	zat_i_lev(:) = 0._wp
[825]	500
[921]	501	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	502	DO ji = 1, nbpac
	503	zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl)
	504	END DO
	505	END DO
[825]	506
[4333]	507	IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
[921]	508	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	509	DO ji = 1, nbpac
[2715]	510	zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
[4333]	511	zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem
	512	zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 )
[2715]	513	END DO
	514	END DO
[4333]	515	IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
[825]	516
[921]	517	!---------------------------------
	518	! Heat content - bottom accretion
	519	!---------------------------------
	520	jm = 1
	521	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	522	DO ji = 1, nbpac
[4161]	523	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise
[2715]	524	zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
	525	zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb &
	526	& + zindb * zdh_frazb(ji) ! frazil ice may coalesce
[3625]	527	zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice
[2715]	528	END DO
	529	END DO
[825]	530
[921]	531	! old layers thicknesses and enthalpies
	532	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	533	DO jk = 1, nlay_i
	534	DO ji = 1, nbpac
[4161]	535	zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i )
[2715]	536	zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
	537	END DO
	538	END DO
	539	END DO
	540	!!gm ??? why the previous do loop if ocerwriten by the following one ?
[921]	541	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	542	DO ji = 1, nbpac
	543	zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl)
[2715]	544	zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl)
[921]	545	END DO ! ji
	546	END DO ! jl
	547
	548	! Redistributing energy on the new grid
[2715]	549	ze_i_ac(:,:,:) = 0._wp
[921]	550	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	551	DO jk = 1, nlay_i
	552	DO layer = 1, nlay_i + 1
	553	DO ji = 1, nbpac
[4161]	554	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )
[921]	555	! Redistributing energy on the new grid
[4161]	556	zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) &
	557	& - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) &
	558	& /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
[2715]	559	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)
[921]	560	END DO ! ji
	561	END DO ! layer
	562	END DO ! jk
	563	END DO ! jl
[825]	564
[921]	565	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
	566	DO jk = 1, nlay_i
	567	DO ji = 1, nbpac
[4161]	568	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )
[2715]	569	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) &
[4161]	570	& / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
[921]	571	END DO
	572	END DO
	573	END DO
[825]	574
[921]	575	!------------
	576	! Update age
	577	!------------
	578	DO jl = 1, jpl
	579	DO ji = 1, nbpac
[4161]	580	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
[2715]	581	zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
	582	END DO
	583	END DO
[825]	584
[921]	585	!-----------------
	586	! Update salinity
	587	!-----------------
[4161]	588	!clem IF( num_sal == 2 ) THEN
[921]	589	DO jl = 1, jpl
	590	DO ji = 1, nbpac
[4161]	591	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
[2715]	592	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
[4161]	593	zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
[2715]	594	END DO
	595	END DO
[4161]	596	!clem ENDIF
[825]	597
[4161]	598	!--------------------------------
	599	! Update mass/salt fluxes (clem)
	600	!--------------------------------
	601	DO jl = 1, jpl
	602	DO ji = 1, nbpac
	603	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
	604	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
	605	rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
	606	sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
	607	END DO
	608	END DO
	609
[921]	610	!------------------------------------------------------------------------------!
	611	! 8) Change 2D vectors to 1D vectors
	612	!------------------------------------------------------------------------------!
	613	DO jl = 1, jpl
[2715]	614	CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
	615	CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
	616	CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
[4161]	617	!clem IF ( num_sal == 2 ) &
[2715]	618	CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
[921]	619	DO jk = 1, nlay_i
[2715]	620	CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
	621	END DO
	622	END DO
[3625]	623	CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
	624	CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
[2715]	625	!
[921]	626	ENDIF ! nbpac > 0
[825]	627
[921]	628	!------------------------------------------------------------------------------!
	629	! 9) Change units for e_i
	630	!------------------------------------------------------------------------------!
[825]	631	DO jl = 1, jpl
[2715]	632	DO jk = 1, nlay_i ! heat content in 10^9 Joules
[4161]	633	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac
[825]	634	END DO
	635	END DO
	636
[921]	637	!------------------------------------------------------------------------------\|
	638	! 10) Conservation check and changes in each ice category
	639	!------------------------------------------------------------------------------\|
[2715]	640	IF( con_i ) THEN
[921]	641	CALL lim_column_sum (jpl, v_i, vt_i_final)
	642	fieldid = 'v_i, limthd_lac'
	643	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
[2715]	644	!
[921]	645	CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
	646	fieldid = 'e_i, limthd_lac'
	647	CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)
[2715]	648	!
[921]	649	CALL lim_column_sum (jpl, v_s, vt_s_final)
	650	fieldid = 'v_s, limthd_lac'
	651	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
[2715]	652	!
[921]	653	! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
	654	! fieldid = 'e_s, limthd_lac'
	655	! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
	656	IF( ln_nicep ) THEN
[4333]	657	DO ji = mi0(jiindx), mi1(jiindx)
	658	DO jj = mj0(jjindx), mj1(jjindx)
	659	WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
	660	WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
	661	WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
	662	WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
	663	END DO
	664	END DO
[921]	665	ENDIF
[2715]	666	!
[834]	667	ENDIF
[2715]	668	!
[3294]	669	CALL wrk_dealloc( jpij, zcatac ) ! integer
	670	CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
	671	CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
	672	CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
	673	CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
	674	CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
	675	CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
[2715]	676	!
[825]	677	END SUBROUTINE lim_thd_lac
	678
	679	#else
[2715]	680	!!----------------------------------------------------------------------
	681	!! Default option NO LIM3 sea-ice model
	682	!!----------------------------------------------------------------------
[825]	683	CONTAINS
	684	SUBROUTINE lim_thd_lac ! Empty routine
	685	END SUBROUTINE lim_thd_lac
	686	#endif
[2715]	687
	688	!!======================================================================
[825]	689	END MODULE limthd_lac

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

Download in other formats: