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.

bdyice.F90 in NEMO/trunk/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/BDY/bdyice.F90 @ 10455

Last change on this file since 10455 was 10425, checked in by smasson, 5 years ago
trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk
Property svn:keywords set to `Id`
File size: 18.2 KB

Rev	Line
[8586]	1	MODULE bdyice
	2	!!======================================================================
	3	!! * MODULE bdyice *
[9656]	4	!! Unstructured Open Boundary Cond. : Open boundary conditions for sea-ice (SI3)
[8586]	5	!!======================================================================
	6	!! History : 3.3 ! 2010-09 (D. Storkey) Original code
[9656]	7	!! 3.4 ! 2012-01 (C. Rousset) add new sea ice model
	8	!! 4.0 ! 2018 (C. Rousset) SI3 compatibility
[8586]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8586]	11	!!----------------------------------------------------------------------
[9656]	12	!! 'key_si3' SI3 sea ice model
[8586]	13	!!----------------------------------------------------------------------
	14	!! bdy_ice : Application of open boundaries to ice
	15	!! bdy_ice_frs : Application of Flow Relaxation Scheme
	16	!!----------------------------------------------------------------------
	17	USE oce ! ocean dynamics and tracers variables
[8637]	18	USE ice ! sea-ice: variables
	19	USE icevar ! sea-ice: operations
[9880]	20	USE icecor ! sea-ice: corrections
[8637]	21	USE icectl ! sea-ice: control prints
[8586]	22	USE phycst ! physical constant
	23	USE eosbn2 ! equation of state
	24	USE par_oce ! ocean parameters
	25	USE dom_oce ! ocean space and time domain variables
	26	USE sbc_oce ! Surface boundary condition: ocean fields
	27	USE bdy_oce ! ocean open boundary conditions
	28	!
	29	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	30	USE in_out_manager ! write to numout file
	31	USE lib_mpp ! distributed memory computing
	32	USE lib_fortran ! to use key_nosignedzero
	33	USE timing ! Timing
	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
	38	PUBLIC bdy_ice ! routine called in sbcmod
[9656]	39	PUBLIC bdy_ice_dyn ! routine called in icedyn_rhg_evp
[8586]	40
	41	!!----------------------------------------------------------------------
[9598]	42	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	43	!! $Id$
[10068]	44	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	45	!!----------------------------------------------------------------------
	46	CONTAINS
	47
	48	SUBROUTINE bdy_ice( kt )
	49	!!----------------------------------------------------------------------
	50	!! * SUBROUTINE bdy_ice *
	51	!!
[9890]	52	!! ** Purpose : Apply open boundary conditions for sea ice
[8586]	53	!!
	54	!!----------------------------------------------------------------------
	55	INTEGER, INTENT(in) :: kt ! Main time step counter
	56	!
[9890]	57	INTEGER :: jbdy ! BDY set index
[8586]	58	!!----------------------------------------------------------------------
	59	!
[9905]	60	IF( ln_timing ) CALL timing_start('bdy_ice_thd')
[9124]	61	!
[8586]	62	CALL ice_var_glo2eqv
	63	!
[9890]	64	DO jbdy = 1, nb_bdy
[8586]	65	!
[9890]	66	SELECT CASE( cn_ice(jbdy) )
[9124]	67	CASE('none') ; CYCLE
[9890]	68	CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy )
[8586]	69	CASE DEFAULT
	70	CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' )
	71	END SELECT
	72	!
	73	END DO
	74	!
[9880]	75	CALL ice_cor( kt , 0 ) ! -- In case categories are out of bounds, do a remapping
[9885]	76	! ! i.e. inputs have not the same ice thickness distribution (set by rn_himean)
	77	! ! than the regional simulation
[9124]	78	CALL ice_var_agg(1)
	79	!
[8586]	80	IF( ln_icectl ) CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
[9905]	81	IF( ln_timing ) CALL timing_stop('bdy_ice_thd')
[8586]	82	!
	83	END SUBROUTINE bdy_ice
	84
	85
[9890]	86	SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy )
[8586]	87	!!------------------------------------------------------------------------------
	88	!! * SUBROUTINE bdy_ice_frs *
	89	!!
[9890]	90	!! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields
[8586]	91	!!
	92	!! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
	93	!! dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
	94	!!------------------------------------------------------------------------------
	95	TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices
	96	TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data
	97	INTEGER, INTENT(in) :: kt ! main time-step counter
[9890]	98	INTEGER, INTENT(in) :: jbdy ! BDY set index
[8586]	99	!
	100	INTEGER :: jpbound ! 0 = incoming ice
	101	! ! 1 = outgoing ice
[9890]	102	INTEGER :: i_bdy, jgrd ! dummy loop indices
	103	INTEGER :: ji, jj, jk, jl, ib, jb
[8586]	104	REAL(wp) :: zwgt, zwgt1 ! local scalar
	105	REAL(wp) :: ztmelts, zdh
	106	!!------------------------------------------------------------------------------
	107	!
	108	jgrd = 1 ! Everything is at T-points here
	109	!
	110	DO jl = 1, jpl
[9890]	111	DO i_bdy = 1, idx%nblenrim(jgrd)
	112	ji = idx%nbi(i_bdy,jgrd)
	113	jj = idx%nbj(i_bdy,jgrd)
	114	zwgt = idx%nbw(i_bdy,jgrd)
	115	zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd)
	116	a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Leads fraction
	117	h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth
	118	h_s(ji,jj,jl) = ( h_s(ji,jj,jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth
[8586]	119
	120	! -----------------
	121	! Pathological case
	122	! -----------------
	123	! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to
	124	! very large transformation from snow to ice (see icethd_dh.F90)
	125
	126	! Then, a) transfer the snow excess into the ice (different from icethd_dh)
[9935]	127	zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )
[8586]	128	! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
[9935]	129	!zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )
[8586]	130
	131	! recompute h_i, h_s
	132	h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
[9935]	133	h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )
[8586]	134
	135	ENDDO
	136	ENDDO
[10425]	137	CALL lbc_bdy_lnk( 'bdyice', a_i(:,:,:), 'T', 1., jbdy )
	138	CALL lbc_bdy_lnk( 'bdyice', h_i(:,:,:), 'T', 1., jbdy )
	139	CALL lbc_bdy_lnk( 'bdyice', h_s(:,:,:), 'T', 1., jbdy )
[8586]	140
	141	DO jl = 1, jpl
[9890]	142	DO i_bdy = 1, idx%nblenrim(jgrd)
	143	ji = idx%nbi(i_bdy,jgrd)
	144	jj = idx%nbj(i_bdy,jgrd)
[8586]	145
	146	! condition on ice thickness depends on the ice velocity
	147	! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
[9890]	148	jpbound = 0 ; ib = ji ; jb = jj
[8586]	149	!
[9890]	150	IF( u_ice(ji+1,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1 ; jb = jj
	151	IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji+1,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1 ; jb = jj
	152	IF( v_ice(ji ,jj+1) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; ib = ji ; jb = jj+1
	153	IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj+1,1) == 0. ) jpbound = 1 ; ib = ji ; jb = jj-1
[8586]	154	!
[9890]	155	IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions
	156	! ! do not make state variables dependent on velocity
[8586]	157	!
[9890]	158	IF( a_i(ib,jb,jl) > 0._wp ) THEN ! there is ice at the boundary
[8586]	159	!
[9890]	160	a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration
	161	h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice
	162	h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw
[8586]	163	!
[9888]	164	SELECT CASE( jpbound )
	165	!
	166	CASE( 0 ) ! velocity is inward
	167	!
[9890]	168	oa_i(ji,jj, jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age
	169	a_ip(ji,jj, jl) = 0._wp ! pond concentration
	170	v_ip(ji,jj, jl) = 0._wp ! pond volume
	171	t_su(ji,jj, jl) = rn_ice_tem(jbdy) ! temperature surface
	172	t_s (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature snw
	173	t_i (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature ice
	174	s_i (ji,jj, jl) = rn_ice_sal(jbdy) ! salinity
	175	sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy) ! salinity profile
[9888]	176	!
	177	CASE( 1 ) ! velocity is outward
	178	!
[9890]	179	oa_i(ji,jj, jl) = oa_i(ib,jb, jl) ! age
	180	a_ip(ji,jj, jl) = a_ip(ib,jb, jl) ! pond concentration
	181	v_ip(ji,jj, jl) = v_ip(ib,jb, jl) ! pond volume
	182	t_su(ji,jj, jl) = t_su(ib,jb, jl) ! temperature surface
	183	t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw
	184	t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice
	185	s_i (ji,jj, jl) = s_i (ib,jb, jl) ! salinity
	186	sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile
[9888]	187	!
	188	END SELECT
[9885]	189	!
[9888]	190	IF( nn_icesal == 1 ) THEN ! if constant salinity
	191	s_i (ji,jj ,jl) = rn_icesal
	192	sz_i(ji,jj,:,jl) = rn_icesal
	193	ENDIF
[8586]	194	!
[9888]	195	! global fields
	196	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) ! volume ice
	197	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) ! volume snw
	198	sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
[8586]	199	DO jk = 1, nlay_s
[9935]	200	e_s(ji,jj,jk,jl) = rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) ! enthalpy in J/m3
[9888]	201	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s ! enthalpy in J/m2
	202	END DO
[8586]	203	DO jk = 1, nlay_i
[9935]	204	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) ! Melting temperature in C
[9888]	205	t_i(ji,jj,jk,jl) = MIN( t_i(ji,jj,jk,jl), ztmelts + rt0 ) ! Force t_i to be lower than melting point => likely conservation issue
	206	!
[9935]	207	e_i(ji,jj,jk,jl) = rhoi * ( rcpi * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) ) & ! enthalpy in J/m3
	208	& + rLfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) ) &
	209	& - rcp * ztmelts )
[9888]	210	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i ! enthalpy in J/m2
[8586]	211	END DO
	212	!
[9888]	213	ELSE ! no ice at the boundary
[9885]	214	!
[9888]	215	a_i (ji,jj, jl) = 0._wp
	216	h_i (ji,jj, jl) = 0._wp
	217	h_s (ji,jj, jl) = 0._wp
	218	oa_i(ji,jj, jl) = 0._wp
	219	a_ip(ji,jj, jl) = 0._wp
	220	v_ip(ji,jj, jl) = 0._wp
	221	t_su(ji,jj, jl) = rt0
	222	t_s (ji,jj,:,jl) = rt0
	223	t_i (ji,jj,:,jl) = rt0
	224
	225	IF( nn_icesal == 1 ) THEN ! if constant salinity
	226	s_i (ji,jj ,jl) = rn_icesal
	227	sz_i(ji,jj,:,jl) = rn_icesal
	228	ELSE ! if variable salinity
	229	s_i (ji,jj,jl) = rn_simin
	230	sz_i(ji,jj,:,jl) = rn_simin
	231	ENDIF
	232	!
	233	! global fields
	234	v_i (ji,jj, jl) = 0._wp
	235	v_s (ji,jj, jl) = 0._wp
	236	sv_i(ji,jj, jl) = 0._wp
	237	e_s (ji,jj,:,jl) = 0._wp
	238	e_i (ji,jj,:,jl) = 0._wp
	239
[8586]	240	ENDIF
[9888]	241
[8586]	242	END DO
	243	!
[9888]	244	END DO ! jl
[9890]	245
[10425]	246	CALL lbc_bdy_lnk( 'bdyice', a_i (:,:,:) , 'T', 1., jbdy )
	247	CALL lbc_bdy_lnk( 'bdyice', h_i (:,:,:) , 'T', 1., jbdy )
	248	CALL lbc_bdy_lnk( 'bdyice', h_s (:,:,:) , 'T', 1., jbdy )
	249	CALL lbc_bdy_lnk( 'bdyice', oa_i(:,:,:) , 'T', 1., jbdy )
	250	CALL lbc_bdy_lnk( 'bdyice', a_ip(:,:,:) , 'T', 1., jbdy )
	251	CALL lbc_bdy_lnk( 'bdyice', v_ip(:,:,:) , 'T', 1., jbdy )
	252	CALL lbc_bdy_lnk( 'bdyice', s_i (:,:,:) , 'T', 1., jbdy )
	253	CALL lbc_bdy_lnk( 'bdyice', t_su(:,:,:) , 'T', 1., jbdy )
	254	CALL lbc_bdy_lnk( 'bdyice', v_i (:,:,:) , 'T', 1., jbdy )
	255	CALL lbc_bdy_lnk( 'bdyice', v_s (:,:,:) , 'T', 1., jbdy )
	256	CALL lbc_bdy_lnk( 'bdyice', sv_i(:,:,:) , 'T', 1., jbdy )
	257	CALL lbc_bdy_lnk( 'bdyice', t_s (:,:,:,:), 'T', 1., jbdy )
	258	CALL lbc_bdy_lnk( 'bdyice', e_s (:,:,:,:), 'T', 1., jbdy )
	259	CALL lbc_bdy_lnk( 'bdyice', t_i (:,:,:,:), 'T', 1., jbdy )
	260	CALL lbc_bdy_lnk( 'bdyice', e_i (:,:,:,:), 'T', 1., jbdy )
[8586]	261	!
	262	END SUBROUTINE bdy_ice_frs
	263
	264
	265	SUBROUTINE bdy_ice_dyn( cd_type )
	266	!!------------------------------------------------------------------------------
	267	!! * SUBROUTINE bdy_ice_dyn *
	268	!!
[9890]	269	!! ** Purpose : Apply dynamics boundary conditions for sea-ice.
[8586]	270	!!
[9890]	271	!! ** Method : if this adjacent grid point is not ice free, then u_ice and v_ice take its value
	272	!! if is ice free, then u_ice and v_ice are unchanged by BDY
	273	!! they keep values calculated in rheology
	274	!!
[8586]	275	!!------------------------------------------------------------------------------
	276	CHARACTER(len=1), INTENT(in) :: cd_type ! nature of velocity grid-points
	277	!
[9890]	278	INTEGER :: i_bdy, jgrd ! dummy loop indices
	279	INTEGER :: ji, jj ! local scalar
	280	INTEGER :: jbdy ! BDY set index
[8586]	281	REAL(wp) :: zmsk1, zmsk2, zflag
	282	!!------------------------------------------------------------------------------
[9905]	283	IF( ln_timing ) CALL timing_start('bdy_ice_dyn')
[8586]	284	!
[9890]	285	DO jbdy=1, nb_bdy
[8586]	286	!
[9890]	287	SELECT CASE( cn_ice(jbdy) )
[8586]	288	!
	289	CASE('none')
	290	CYCLE
	291	!
	292	CASE('frs')
	293	!
[9890]	294	IF( nn_ice_dta(jbdy) == 0 ) CYCLE ! case ice boundaries = initial conditions
	295	! ! do not change ice velocity (it is only computed by rheology)
[8586]	296	SELECT CASE ( cd_type )
	297	!
	298	CASE ( 'U' )
	299	jgrd = 2 ! u velocity
[9890]	300	DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd)
	301	ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
	302	jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
	303	zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd)
[8586]	304	!
	305	IF ( ABS( zflag ) == 1. ) THEN ! eastern and western boundaries
	306	! one of the two zmsk is always 0 (because of zflag)
	307	zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
	308	zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
	309	!
[9888]	310	! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice)
[8586]	311	u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
	312	& u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
[9888]	313	& u_ice(ji ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
[8586]	314	ELSE ! everywhere else
	315	u_ice(ji,jj) = 0._wp
	316	ENDIF
	317	!
	318	END DO
[10425]	319	CALL lbc_bdy_lnk( 'bdyice', u_ice(:,:), 'U', -1., jbdy )
[8586]	320	!
	321	CASE ( 'V' )
	322	jgrd = 3 ! v velocity
[9890]	323	DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd)
	324	ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
	325	jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
	326	zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd)
[8586]	327	!
	328	IF ( ABS( zflag ) == 1. ) THEN ! northern and southern boundaries
	329	! one of the two zmsk is always 0 (because of zflag)
	330	zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
	331	zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
	332	!
[9888]	333	! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice)
[8586]	334	v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
	335	& v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
[9888]	336	& v_ice(ji,jj ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
[8586]	337	ELSE ! everywhere else
	338	v_ice(ji,jj) = 0._wp
	339	ENDIF
	340	!
	341	END DO
[10425]	342	CALL lbc_bdy_lnk( 'bdyice', v_ice(:,:), 'V', -1., jbdy )
[8586]	343	!
	344	END SELECT
	345	!
	346	CASE DEFAULT
	347	CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' )
	348	END SELECT
	349	!
	350	END DO
	351	!
[9905]	352	IF( ln_timing ) CALL timing_stop('bdy_ice_dyn')
	353	!
[8586]	354	END SUBROUTINE bdy_ice_dyn
	355
	356	#else
	357	!!---------------------------------------------------------------------------------
	358	!! Default option Empty module
	359	!!---------------------------------------------------------------------------------
	360	CONTAINS
	361	SUBROUTINE bdy_ice( kt ) ! Empty routine
[9927]	362	IMPLICIT NONE
	363	INTEGER, INTENT( in ) :: kt
[8586]	364	WRITE(,) 'bdy_ice: You should not have seen this print! error?', kt
	365	END SUBROUTINE bdy_ice
	366	#endif
	367
	368	!!=================================================================================
	369	END MODULE bdyice

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

Download in other formats: