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.

icethd_do.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icethd_do.F90 @ 13601

Last change on this file since 13601 was 13601, checked in by clem, 4 years ago
trunk: rewrite heat budget of sea ice to make it perfectly conservative by construction. Also, activating ln_icediachk now gives an ascii file (icedrift_diagnostics.ascii) containing mass, salt and heat global conservation issues (if any). In addition, 2D drift files can be outputed (icedrift_mass…) and field_def is changed accordingly. Note that advection schemes are not yet commited since there is still a restartability issue that I do not understand.
Property svn:keywords set to `Id`
File size: 23.2 KB

Rev	Line
[8586]	1	MODULE icethd_do
	2	!!======================================================================
	3	!! * MODULE icethd_do *
	4	!! sea-ice: sea ice growth in the leads (open water)
	5	!!======================================================================
[9604]	6	!! History : ! 2005-12 (M. Vancoppenolle) Original code
	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	8	!!----------------------------------------------------------------------
[9570]	9	#if defined key_si3
[8586]	10	!!----------------------------------------------------------------------
[9570]	11	!! 'key_si3' SI3 sea-ice model
[8586]	12	!!----------------------------------------------------------------------
	13	!! ice_thd_do : ice growth in open water (=lateral accretion of ice)
	14	!! ice_thd_do_init : initialization
	15	!!----------------------------------------------------------------------
	16	USE dom_oce ! ocean space and time domain
	17	USE phycst ! physical constants
	18	USE sbc_oce , ONLY : sss_m
	19	USE sbc_ice , ONLY : utau_ice, vtau_ice
	20	USE ice1D ! sea-ice: thermodynamics variables
	21	USE ice ! sea-ice: variables
	22	USE icetab ! sea-ice: 2D <==> 1D
	23	USE icectl ! sea-ice: conservation
	24	USE icethd_ent ! sea-ice: thermodynamics, enthalpy
	25	USE icevar ! sea-ice: operations
	26	USE icethd_sal ! sea-ice: salinity profiles
	27	!
	28	USE in_out_manager ! I/O manager
	29	USE lib_mpp ! MPP library
	30	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	31	USE lbclnk ! lateral boundary conditions (or mpp links)
	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC ice_thd_do ! called by ice_thd
	37	PUBLIC ice_thd_do_init ! called by ice_stp
	38
[9169]	39	! !! namelist (namthd_do)
	40	REAL(wp) :: rn_hinew ! thickness for new ice formation (m)
	41	LOGICAL :: ln_frazil ! use of frazil ice collection as function of wind (T) or not (F)
	42	REAL(wp) :: rn_maxfraz ! maximum portion of frazil ice collecting at the ice bottom
	43	REAL(wp) :: rn_vfraz ! threshold drift speed for collection of bottom frazil ice
	44	REAL(wp) :: rn_Cfraz ! squeezing coefficient for collection of bottom frazil ice
[8586]	45
[12377]	46	!! * Substitutions
	47	# include "do_loop_substitute.h90"
[8586]	48	!!----------------------------------------------------------------------
[9598]	49	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	50	!! $Id$
[10068]	51	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	52	!!----------------------------------------------------------------------
	53	CONTAINS
	54
	55	SUBROUTINE ice_thd_do
	56	!!-------------------------------------------------------------------
	57	!! * ROUTINE ice_thd_do *
	58	!!
	59	!! ** Purpose : Computation of the evolution of the ice thickness and
[9604]	60	!! concentration as a function of the heat balance in the leads
[8586]	61	!!
[9604]	62	!! ** Method : Ice is formed in the open water when ocean looses heat
	63	!! (heat budget of open water is negative) following
	64	!!
	65	!! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
	66	!! where - h0 is the thickness of ice created in the lead
	67	!! - a is a minimum fraction for leads
	68	!! - F is a monotonic non-increasing function defined as:
[8586]	69	!! F(X)=( 1 - Xexld )(1.0/exld)
[9604]	70	!! - exld is the exponent closure rate (=2 default val.)
[8586]	71	!!
	72	!! ** Action : - Adjustment of snow and ice thicknesses and heat
	73	!! content in brine pockets
	74	!! - Updating ice internal temperature
	75	!! - Computation of variation of ice volume and mass
	76	!! - Computation of a_i after lateral accretion and
	77	!! update h_s_1d, h_i_1d
	78	!!------------------------------------------------------------------------
[9169]	79	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[9604]	80	INTEGER :: iter ! - -
	81	REAL(wp) :: ztmelts, zfrazb, zweight, zde ! local scalars
[8586]	82	REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf ! - -
	83	REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - -
[9169]	84	!
[8586]	85	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
	86	REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg)
	87	REAL(wp) :: zEw ! seawater specific enthalpy (J/kg)
	88	REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean)
[9169]	89	!
[8586]	90	REAL(wp) :: zv_newfra
[9169]	91	!
[8586]	92	INTEGER , DIMENSION(jpij) :: jcat ! indexes of categories where new ice grows
	93	REAL(wp), DIMENSION(jpij) :: zswinew ! switch for new ice or not
[9169]	94	!
[8586]	95	REAL(wp), DIMENSION(jpij) :: zv_newice ! volume of accreted ice
	96	REAL(wp), DIMENSION(jpij) :: za_newice ! fractional area of accreted ice
	97	REAL(wp), DIMENSION(jpij) :: zh_newice ! thickness of accreted ice
	98	REAL(wp), DIMENSION(jpij) :: ze_newice ! heat content of accreted ice
	99	REAL(wp), DIMENSION(jpij) :: zs_newice ! salinity of accreted ice
	100	REAL(wp), DIMENSION(jpij) :: zo_newice ! age of accreted ice
	101	REAL(wp), DIMENSION(jpij) :: zdv_res ! residual volume in case of excessive heat budget
	102	REAL(wp), DIMENSION(jpij) :: zda_res ! residual area in case of excessive heat budget
	103	REAL(wp), DIMENSION(jpij) :: zv_frazb ! accretion of frazil ice at the ice bottom
	104	REAL(wp), DIMENSION(jpij) :: zvrel_1d ! relative ice / frazil velocity (1D vector)
[9169]	105	!
[9604]	106	REAL(wp), DIMENSION(jpij,jpl) :: zv_b ! old volume of ice in category jl
	107	REAL(wp), DIMENSION(jpij,jpl) :: za_b ! old area of ice in category jl
[9169]	108	!
[8586]	109	REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d !: 1-D version of e_i
[9169]	110	!
[9604]	111	REAL(wp), DIMENSION(jpi,jpj) :: zvrel ! relative ice / frazil velocity
[9169]	112	!
[9604]	113	REAL(wp) :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used)
[8586]	114	!!-----------------------------------------------------------------------!
	115
[9169]	116	IF( ln_icediachk ) CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
[11536]	117	IF( ln_icediachk ) CALL ice_cons2D ( 0, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft )
[8586]	118
[10994]	119	at_i(:,:) = SUM( a_i, dim=3 )
[8586]	120	!------------------------------------------------------------------------------!
[9604]	121	! 1) Collection thickness of ice formed in leads and polynyas
[8586]	122	!------------------------------------------------------------------------------!
	123	! ht_i_new is the thickness of new ice formed in open water
[9604]	124	! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
[8586]	125	! Frazil ice forms in open water, is transported by wind
	126	! accumulates at the edge of the consolidated ice edge
	127	! where it forms aggregates of a specific thickness called
	128	! collection thickness.
	129
	130	zvrel(:,:) = 0._wp
	131
	132	! Default new ice thickness
[13601]	133	WHERE( qlead(:,:) < 0._wp ) ! cooling
	134	ht_i_new(:,:) = rn_hinew
	135	ELSEWHERE
	136	ht_i_new(:,:) = 0._wp
[8586]	137	END WHERE
	138
	139	IF( ln_frazil ) THEN
[9169]	140	!
[8586]	141	ht_i_new(:,:) = 0._wp
[9169]	142	!
[9604]	143	! Physical constants
	144	zhicrit = 0.04 ! frazil ice thickness
[12489]	145	ztwogp = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) ) ! reduced grav
[9604]	146	zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)
[8586]	147	zgamafr = 0.03
[9169]	148	!
[13295]	149	DO_2D( 0, 0, 0, 0 )
[13601]	150	IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
[12377]	151	! -- Wind stress -- !
	152	ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) &
	153	& + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp
	154	ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) &
	155	& + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp
	156	! Square root of wind stress
	157	ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
[8586]	158
[12377]	159	! -- Frazil ice velocity -- !
	160	rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
	161	zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
	162	zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
[8586]	163
[12377]	164	! -- Pack ice velocity -- !
	165	zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
	166	zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
[8586]	167
[12377]	168	! -- Relative frazil/pack ice velocity -- !
	169	rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
	170	zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) &
	171	& + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
	172	zvrel(ji,jj) = SQRT( zvrel2 )
[8586]	173
[12377]	174	! -- new ice thickness (iterative loop) -- !
	175	ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1 ) &
	176	& / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2
[8586]	177
[12377]	178	iter = 1
	179	DO WHILE ( iter < 20 )
	180	zf = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) - &
	181	& ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
	182	zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
[8586]	183
[12377]	184	ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
	185	iter = iter + 1
	186	END DO
[9169]	187	!
[12377]	188	! bound ht_i_new (though I don't see why it should be necessary)
	189	ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
	190	!
	191	ENDIF
	192	!
	193	END_2D
[8586]	194	!
[13226]	195	CALL lbc_lnk_multi( 'icethd_do', zvrel, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp )
[8586]	196
[9604]	197	ENDIF
[8586]	198
	199	!------------------------------------------------------------------------------!
[9604]	200	! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
[8586]	201	!------------------------------------------------------------------------------!
[13601]	202	! it occurs if cooling
[8586]	203
[9604]	204	! Identify grid points where new ice forms
[9169]	205	npti = 0 ; nptidx(:) = 0
[13295]	206	DO_2D( 1, 1, 1, 1 )
[13601]	207	IF ( qlead(ji,jj) < 0._wp ) THEN
[12377]	208	npti = npti + 1
	209	nptidx( npti ) = (jj - 1) * jpi + ji
	210	ENDIF
	211	END_2D
[8586]	212
	213	! Move from 2-D to 1-D vectors
	214	IF ( npti > 0 ) THEN
	215
	216	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , at_i )
	217	CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
	218	CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
	219	CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
	220	DO jl = 1, jpl
	221	DO jk = 1, nlay_i
	222	CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
	223	END DO
	224	END DO
[9604]	225	CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead )
	226	CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo )
	227	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw )
	228	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw )
	229	CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new )
	230	CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d (1:npti) , zvrel )
[8586]	231
[9604]	232	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd )
	233	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw )
	234	CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
	235	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m )
[8586]	236
[9604]	237	! Convert units for ice internal energy
[8586]	238	DO jl = 1, jpl
	239	DO jk = 1, nlay_i
	240	WHERE( v_i_2d(1:npti,jl) > 0._wp )
	241	ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
	242	ELSEWHERE
	243	ze_i_2d(1:npti,jk,jl) = 0._wp
	244	END WHERE
	245	END DO
	246	END DO
	247
	248	! Keep old ice areas and volume in memory
	249	zv_b(1:npti,:) = v_i_2d(1:npti,:)
	250	za_b(1:npti,:) = a_i_2d(1:npti,:)
	251
[9604]	252	! --- Salinity of new ice --- !
[8586]	253	SELECT CASE ( nn_icesal )
	254	CASE ( 1 ) ! Sice = constant
	255	zs_newice(1:npti) = rn_icesal
	256	CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
	257	DO ji = 1, npti
	258	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
	259	END DO
	260	CASE ( 3 ) ! Sice = F(z) [multiyear ice]
	261	zs_newice(1:npti) = 2.3
	262	END SELECT
	263
[9604]	264	! --- Heat content of new ice --- !
[8586]	265	! We assume that new ice is formed at the seawater freezing point
	266	DO ji = 1, npti
[9935]	267	ztmelts = - rTmlt * zs_newice(ji) ! Melting point (C)
	268	ze_newice(ji) = rhoi * ( rcpi * ( ztmelts - ( t_bo_1d(ji) - rt0 ) ) &
	269	& + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) &
	270	& - rcp * ztmelts )
[8586]	271	END DO
	272
[9604]	273	! --- Age of new ice --- !
[8586]	274	zo_newice(1:npti) = 0._wp
	275
[9604]	276	! --- Volume of new ice --- !
[8586]	277	DO ji = 1, npti
	278
[9935]	279	zEi = - ze_newice(ji) * r1_rhoi ! specific enthalpy of forming ice [J/kg]
[8586]	280
	281	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg]
	282	! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied)
	283
	284	zdE = zEi - zEw ! specific enthalpy difference [J/kg]
	285
	286	zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0)
	287	! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point
[9935]	288	zv_newice(ji) = - zfmdt * r1_rhoi
[8586]	289
	290	zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux
	291
	292	! Contribution to heat flux to the ocean [W.m-2], >0
[12489]	293	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_Dt_ice
[8586]	294	! Total heat flux used in this process [W.m-2]
[12489]	295	hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice
[8586]	296	! mass flux
[12489]	297	wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_Dt_ice
[8586]	298	! salt flux
[12489]	299	sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice
[8586]	300	END DO
	301
	302	zv_frazb(1:npti) = 0._wp
	303	IF( ln_frazil ) THEN
	304	! A fraction zfrazb of frazil ice is accreted at the ice bottom
	305	DO ji = 1, npti
	306	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
	307	zfrazb = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
	308	zv_frazb(ji) = zfrazb * zv_newice(ji)
	309	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
	310	END DO
	311	END IF
	312
[9604]	313	! --- Area of new ice --- !
[8586]	314	DO ji = 1, npti
	315	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
	316	END DO
	317
	318	!------------------------------------------------------------------------------!
[9604]	319	! 3) Redistribute new ice area and volume into ice categories !
[8586]	320	!------------------------------------------------------------------------------!
	321
[9604]	322	! --- lateral ice growth --- !
[10994]	323	! If lateral ice growth gives an ice concentration > amax, then
[8586]	324	! we keep the excessive volume in memory and attribute it later to bottom accretion
	325	DO ji = 1, npti
[10994]	326	IF ( za_newice(ji) > MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN ! max is for roundoff error
	327	zda_res(ji) = za_newice(ji) - MAX( 0._wp, rn_amax_1d(ji) - at_i_1d(ji) )
[8586]	328	zdv_res(ji) = zda_res (ji) * zh_newice(ji)
[10994]	329	za_newice(ji) = MAX( 0._wp, za_newice(ji) - zda_res (ji) )
	330	zv_newice(ji) = MAX( 0._wp, zv_newice(ji) - zdv_res (ji) )
[8586]	331	ELSE
	332	zda_res(ji) = 0._wp
	333	zdv_res(ji) = 0._wp
	334	ENDIF
	335	END DO
	336
	337	! find which category to fill
	338	DO jl = 1, jpl
	339	DO ji = 1, npti
	340	IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
	341	a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
	342	v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
	343	jcat(ji) = jl
	344	ENDIF
	345	END DO
	346	END DO
	347	at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
	348
	349	! Heat content
	350	DO ji = 1, npti
	351	jl = jcat(ji) ! categroy in which new ice is put
	352	zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice
	353	END DO
	354
	355	DO jk = 1, nlay_i
	356	DO ji = 1, npti
	357	jl = jcat(ji)
	358	rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
	359	ze_i_2d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + &
	360	& ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) ) &
	361	& * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
	362	END DO
	363	END DO
	364
[9604]	365	! --- bottom ice growth + ice enthalpy remapping --- !
[8586]	366	DO jl = 1, jpl
	367
	368	! for remapping
	369	h_i_old (1:npti,0:nlay_i+1) = 0._wp
	370	eh_i_old(1:npti,0:nlay_i+1) = 0._wp
	371	DO jk = 1, nlay_i
	372	DO ji = 1, npti
	373	h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
	374	eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
	375	END DO
	376	END DO
	377
	378	! new volumes including lateral/bottom accretion + residual
	379	DO ji = 1, npti
	380	rswitch = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
	381	zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
	382	a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)
	383	v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
	384	! for remapping
	385	h_i_old (ji,nlay_i+1) = zv_newfra
	386	eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
[9169]	387	END DO
[8586]	388	! --- Ice enthalpy remapping --- !
[13547]	389	CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) )
[9169]	390	END DO
[8586]	391
[9604]	392	! --- Update salinity --- !
[8586]	393	DO jl = 1, jpl
	394	DO ji = 1, npti
	395	sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
	396	END DO
	397	END DO
	398
[9604]	399	! Change units for e_i
[8586]	400	DO jl = 1, jpl
	401	DO jk = 1, nlay_i
	402	ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i
	403	END DO
	404	END DO
[9604]	405
	406	! Move 2D vectors to 1D vectors
[8586]	407	CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
	408	CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
	409	CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
	410	DO jl = 1, jpl
	411	DO jk = 1, nlay_i
	412	CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
	413	END DO
	414	END DO
	415	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
	416	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
	417	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
	418	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
	419	!
	420	ENDIF ! npti > 0
	421	!
	422	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
[11536]	423	IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd_do', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft)
[8586]	424	!
	425	END SUBROUTINE ice_thd_do
	426
[9169]	427
[8586]	428	SUBROUTINE ice_thd_do_init
	429	!!-----------------------------------------------------------------------
	430	!! * ROUTINE ice_thd_do_init *
	431	!!
	432	!! ** Purpose : Physical constants and parameters associated with
	433	!! ice growth in the leads
	434	!!
	435	!! ** Method : Read the namthd_do namelist and check the parameters
	436	!! called at the first timestep (nit000)
	437	!!
	438	!! ** input : Namelist namthd_do
	439	!!-------------------------------------------------------------------
[9169]	440	INTEGER :: ios ! Local integer
[8586]	441	!!
	442	NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
	443	!!-------------------------------------------------------------------
	444	!
	445	READ ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
[11536]	446	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_do in reference namelist' )
[8586]	447	READ ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
[11536]	448	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_do in configuration namelist' )
[9169]	449	IF(lwm) WRITE( numoni, namthd_do )
[8586]	450	!
	451	IF(lwp) THEN ! control print
[9169]	452	WRITE(numout,*)
[8586]	453	WRITE(numout,*) 'ice_thd_do_init: Ice growth in open water'
	454	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	455	WRITE(numout,*) ' Namelist namthd_do:'
	456	WRITE(numout,*) ' ice thickness for lateral accretion rn_hinew = ', rn_hinew
	457	WRITE(numout,*) ' Frazil ice thickness as a function of wind or not ln_frazil = ', ln_frazil
	458	WRITE(numout,*) ' Maximum proportion of frazil ice collecting at bottom rn_maxfraz = ', rn_maxfraz
	459	WRITE(numout,*) ' Threshold relative drift speed for collection of frazil rn_vfraz = ', rn_vfraz
	460	WRITE(numout,*) ' Squeezing coefficient for collection of frazil rn_Cfraz = ', rn_Cfraz
	461	ENDIF
	462	!
	463	IF ( rn_hinew < rn_himin ) CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )
	464	!
	465	END SUBROUTINE ice_thd_do_init
	466
	467	#else
	468	!!----------------------------------------------------------------------
[9570]	469	!! Default option NO SI3 sea-ice model
[8586]	470	!!----------------------------------------------------------------------
	471	#endif
	472
	473	!!======================================================================
	474	END MODULE icethd_do

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

Download in other formats: