Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

bdyice.F90 @ 11083

Last change on this file since 11083 was 11083, checked in by davestorkey, 5 years ago
UKMO/NEMO_4.0_GO6_mixing : Update to be relative to rev 11081 of NEMO_4.0_mirror.
File size: 18.6 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	!!----------------------------------------------------------------------
[11083]	59	! controls
	60	IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing
	61	IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[8586]	62	!
	63	CALL ice_var_glo2eqv
	64	!
[9890]	65	DO jbdy = 1, nb_bdy
[8586]	66	!
[9890]	67	SELECT CASE( cn_ice(jbdy) )
[9124]	68	CASE('none') ; CYCLE
[9890]	69	CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy )
[8586]	70	CASE DEFAULT
	71	CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' )
	72	END SELECT
	73	!
	74	END DO
	75	!
[9880]	76	CALL ice_cor( kt , 0 ) ! -- In case categories are out of bounds, do a remapping
[9885]	77	! ! i.e. inputs have not the same ice thickness distribution (set by rn_himean)
	78	! ! than the regional simulation
[9124]	79	CALL ice_var_agg(1)
	80	!
[11083]	81	! controls
	82	IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
	83	IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints
	84	IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing
[8586]	85	!
	86	END SUBROUTINE bdy_ice
	87
	88
[9890]	89	SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy )
[8586]	90	!!------------------------------------------------------------------------------
	91	!! * SUBROUTINE bdy_ice_frs *
	92	!!
[9890]	93	!! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields
[8586]	94	!!
	95	!! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
	96	!! dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
	97	!!------------------------------------------------------------------------------
	98	TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices
	99	TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data
	100	INTEGER, INTENT(in) :: kt ! main time-step counter
[9890]	101	INTEGER, INTENT(in) :: jbdy ! BDY set index
[8586]	102	!
	103	INTEGER :: jpbound ! 0 = incoming ice
	104	! ! 1 = outgoing ice
[9890]	105	INTEGER :: i_bdy, jgrd ! dummy loop indices
	106	INTEGER :: ji, jj, jk, jl, ib, jb
[8586]	107	REAL(wp) :: zwgt, zwgt1 ! local scalar
	108	REAL(wp) :: ztmelts, zdh
	109	!!------------------------------------------------------------------------------
	110	!
	111	jgrd = 1 ! Everything is at T-points here
	112	!
	113	DO jl = 1, jpl
[9890]	114	DO i_bdy = 1, idx%nblenrim(jgrd)
	115	ji = idx%nbi(i_bdy,jgrd)
	116	jj = idx%nbj(i_bdy,jgrd)
	117	zwgt = idx%nbw(i_bdy,jgrd)
	118	zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd)
	119	a_i(ji,jj,jl) = ( a_i(ji,jj,jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Leads fraction
	120	h_i(ji,jj,jl) = ( h_i(ji,jj,jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth
	121	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]	122
	123	! -----------------
	124	! Pathological case
	125	! -----------------
	126	! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to
	127	! very large transformation from snow to ice (see icethd_dh.F90)
	128
	129	! Then, a) transfer the snow excess into the ice (different from icethd_dh)
[9935]	130	zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )
[8586]	131	! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
[9935]	132	!zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )
[8586]	133
	134	! recompute h_i, h_s
	135	h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
[9935]	136	h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )
[8586]	137
	138	ENDDO
	139	ENDDO
[10425]	140	CALL lbc_bdy_lnk( 'bdyice', a_i(:,:,:), 'T', 1., jbdy )
	141	CALL lbc_bdy_lnk( 'bdyice', h_i(:,:,:), 'T', 1., jbdy )
	142	CALL lbc_bdy_lnk( 'bdyice', h_s(:,:,:), 'T', 1., jbdy )
[8586]	143
	144	DO jl = 1, jpl
[9890]	145	DO i_bdy = 1, idx%nblenrim(jgrd)
	146	ji = idx%nbi(i_bdy,jgrd)
	147	jj = idx%nbj(i_bdy,jgrd)
[8586]	148
	149	! condition on ice thickness depends on the ice velocity
	150	! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
[9890]	151	jpbound = 0 ; ib = ji ; jb = jj
[8586]	152	!
[11083]	153	IF( u_ice(ji ,jj ) < 0. .AND. umask(ji-1,jj ,1) == 0. ) jpbound = 1 ; ib = ji+1
	154	IF( u_ice(ji-1,jj ) > 0. .AND. umask(ji ,jj ,1) == 0. ) jpbound = 1 ; ib = ji-1
	155	IF( v_ice(ji ,jj ) < 0. .AND. vmask(ji ,jj-1,1) == 0. ) jpbound = 1 ; jb = jj+1
	156	IF( v_ice(ji ,jj-1) > 0. .AND. vmask(ji ,jj ,1) == 0. ) jpbound = 1 ; jb = jj-1
[8586]	157	!
[9890]	158	IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions
	159	! ! do not make state variables dependent on velocity
[8586]	160	!
[9890]	161	IF( a_i(ib,jb,jl) > 0._wp ) THEN ! there is ice at the boundary
[8586]	162	!
[9890]	163	a_i(ji,jj,jl) = a_i(ib,jb,jl) ! concentration
	164	h_i(ji,jj,jl) = h_i(ib,jb,jl) ! thickness ice
	165	h_s(ji,jj,jl) = h_s(ib,jb,jl) ! thickness snw
[8586]	166	!
[9888]	167	SELECT CASE( jpbound )
	168	!
	169	CASE( 0 ) ! velocity is inward
	170	!
[9890]	171	oa_i(ji,jj, jl) = rn_ice_age(jbdy) * a_i(ji,jj,jl) ! age
	172	a_ip(ji,jj, jl) = 0._wp ! pond concentration
	173	v_ip(ji,jj, jl) = 0._wp ! pond volume
	174	t_su(ji,jj, jl) = rn_ice_tem(jbdy) ! temperature surface
	175	t_s (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature snw
	176	t_i (ji,jj,:,jl) = rn_ice_tem(jbdy) ! temperature ice
	177	s_i (ji,jj, jl) = rn_ice_sal(jbdy) ! salinity
	178	sz_i(ji,jj,:,jl) = rn_ice_sal(jbdy) ! salinity profile
[9888]	179	!
	180	CASE( 1 ) ! velocity is outward
	181	!
[9890]	182	oa_i(ji,jj, jl) = oa_i(ib,jb, jl) ! age
	183	a_ip(ji,jj, jl) = a_ip(ib,jb, jl) ! pond concentration
	184	v_ip(ji,jj, jl) = v_ip(ib,jb, jl) ! pond volume
	185	t_su(ji,jj, jl) = t_su(ib,jb, jl) ! temperature surface
	186	t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl) ! temperature snw
	187	t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl) ! temperature ice
	188	s_i (ji,jj, jl) = s_i (ib,jb, jl) ! salinity
	189	sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl) ! salinity profile
[9888]	190	!
	191	END SELECT
[9885]	192	!
[9888]	193	IF( nn_icesal == 1 ) THEN ! if constant salinity
	194	s_i (ji,jj ,jl) = rn_icesal
	195	sz_i(ji,jj,:,jl) = rn_icesal
	196	ENDIF
[8586]	197	!
[9888]	198	! global fields
	199	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) ! volume ice
	200	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) ! volume snw
	201	sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
[8586]	202	DO jk = 1, nlay_s
[9935]	203	e_s(ji,jj,jk,jl) = rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) ! enthalpy in J/m3
[9888]	204	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s ! enthalpy in J/m2
	205	END DO
[8586]	206	DO jk = 1, nlay_i
[9935]	207	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) ! Melting temperature in C
[9888]	208	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
	209	!
[9935]	210	e_i(ji,jj,jk,jl) = rhoi * ( rcpi * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) ) & ! enthalpy in J/m3
	211	& + rLfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) ) &
	212	& - rcp * ztmelts )
[9888]	213	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]	214	END DO
	215	!
[9888]	216	ELSE ! no ice at the boundary
[9885]	217	!
[9888]	218	a_i (ji,jj, jl) = 0._wp
	219	h_i (ji,jj, jl) = 0._wp
	220	h_s (ji,jj, jl) = 0._wp
	221	oa_i(ji,jj, jl) = 0._wp
	222	a_ip(ji,jj, jl) = 0._wp
	223	v_ip(ji,jj, jl) = 0._wp
	224	t_su(ji,jj, jl) = rt0
	225	t_s (ji,jj,:,jl) = rt0
	226	t_i (ji,jj,:,jl) = rt0
	227
	228	IF( nn_icesal == 1 ) THEN ! if constant salinity
	229	s_i (ji,jj ,jl) = rn_icesal
	230	sz_i(ji,jj,:,jl) = rn_icesal
	231	ELSE ! if variable salinity
	232	s_i (ji,jj,jl) = rn_simin
	233	sz_i(ji,jj,:,jl) = rn_simin
	234	ENDIF
	235	!
	236	! global fields
	237	v_i (ji,jj, jl) = 0._wp
	238	v_s (ji,jj, jl) = 0._wp
	239	sv_i(ji,jj, jl) = 0._wp
	240	e_s (ji,jj,:,jl) = 0._wp
	241	e_i (ji,jj,:,jl) = 0._wp
	242
[8586]	243	ENDIF
[9888]	244
[8586]	245	END DO
	246	!
[9888]	247	END DO ! jl
[9890]	248
[10425]	249	CALL lbc_bdy_lnk( 'bdyice', a_i (:,:,:) , 'T', 1., jbdy )
	250	CALL lbc_bdy_lnk( 'bdyice', h_i (:,:,:) , 'T', 1., jbdy )
	251	CALL lbc_bdy_lnk( 'bdyice', h_s (:,:,:) , 'T', 1., jbdy )
	252	CALL lbc_bdy_lnk( 'bdyice', oa_i(:,:,:) , 'T', 1., jbdy )
	253	CALL lbc_bdy_lnk( 'bdyice', a_ip(:,:,:) , 'T', 1., jbdy )
	254	CALL lbc_bdy_lnk( 'bdyice', v_ip(:,:,:) , 'T', 1., jbdy )
	255	CALL lbc_bdy_lnk( 'bdyice', s_i (:,:,:) , 'T', 1., jbdy )
	256	CALL lbc_bdy_lnk( 'bdyice', t_su(:,:,:) , 'T', 1., jbdy )
	257	CALL lbc_bdy_lnk( 'bdyice', v_i (:,:,:) , 'T', 1., jbdy )
	258	CALL lbc_bdy_lnk( 'bdyice', v_s (:,:,:) , 'T', 1., jbdy )
	259	CALL lbc_bdy_lnk( 'bdyice', sv_i(:,:,:) , 'T', 1., jbdy )
	260	CALL lbc_bdy_lnk( 'bdyice', t_s (:,:,:,:), 'T', 1., jbdy )
	261	CALL lbc_bdy_lnk( 'bdyice', e_s (:,:,:,:), 'T', 1., jbdy )
	262	CALL lbc_bdy_lnk( 'bdyice', t_i (:,:,:,:), 'T', 1., jbdy )
	263	CALL lbc_bdy_lnk( 'bdyice', e_i (:,:,:,:), 'T', 1., jbdy )
[8586]	264	!
	265	END SUBROUTINE bdy_ice_frs
	266
	267
	268	SUBROUTINE bdy_ice_dyn( cd_type )
	269	!!------------------------------------------------------------------------------
	270	!! * SUBROUTINE bdy_ice_dyn *
	271	!!
[9890]	272	!! ** Purpose : Apply dynamics boundary conditions for sea-ice.
[8586]	273	!!
[9890]	274	!! ** Method : if this adjacent grid point is not ice free, then u_ice and v_ice take its value
	275	!! if is ice free, then u_ice and v_ice are unchanged by BDY
	276	!! they keep values calculated in rheology
	277	!!
[8586]	278	!!------------------------------------------------------------------------------
	279	CHARACTER(len=1), INTENT(in) :: cd_type ! nature of velocity grid-points
	280	!
[9890]	281	INTEGER :: i_bdy, jgrd ! dummy loop indices
	282	INTEGER :: ji, jj ! local scalar
	283	INTEGER :: jbdy ! BDY set index
[8586]	284	REAL(wp) :: zmsk1, zmsk2, zflag
	285	!!------------------------------------------------------------------------------
[9905]	286	IF( ln_timing ) CALL timing_start('bdy_ice_dyn')
[8586]	287	!
[9890]	288	DO jbdy=1, nb_bdy
[8586]	289	!
[9890]	290	SELECT CASE( cn_ice(jbdy) )
[8586]	291	!
	292	CASE('none')
	293	CYCLE
	294	!
	295	CASE('frs')
	296	!
[9890]	297	IF( nn_ice_dta(jbdy) == 0 ) CYCLE ! case ice boundaries = initial conditions
	298	! ! do not change ice velocity (it is only computed by rheology)
[8586]	299	SELECT CASE ( cd_type )
	300	!
	301	CASE ( 'U' )
	302	jgrd = 2 ! u velocity
[9890]	303	DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd)
	304	ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
	305	jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
	306	zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd)
[8586]	307	!
	308	IF ( ABS( zflag ) == 1. ) THEN ! eastern and western boundaries
	309	! one of the two zmsk is always 0 (because of zflag)
	310	zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
[11083]	311	zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) ) ! 0 if no ice
[8586]	312	!
[9888]	313	! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then do not change u_ice)
[8586]	314	u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
	315	& u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
[9888]	316	& u_ice(ji ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
[8586]	317	ELSE ! everywhere else
	318	u_ice(ji,jj) = 0._wp
	319	ENDIF
	320	!
	321	END DO
[10425]	322	CALL lbc_bdy_lnk( 'bdyice', u_ice(:,:), 'U', -1., jbdy )
[8586]	323	!
	324	CASE ( 'V' )
	325	jgrd = 3 ! v velocity
[9890]	326	DO i_bdy = 1, idx_bdy(jbdy)%nblenrim(jgrd)
	327	ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
	328	jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
	329	zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd)
[8586]	330	!
	331	IF ( ABS( zflag ) == 1. ) THEN ! northern and southern boundaries
	332	! one of the two zmsk is always 0 (because of zflag)
	333	zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
[11083]	334	zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj) ) ) ! 0 if no ice
[8586]	335	!
[9888]	336	! v_ice = v_ice of the adjacent grid point except if this grid point is ice-free (then do not change v_ice)
[8586]	337	v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
	338	& v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
[9888]	339	& v_ice(ji,jj ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
[8586]	340	ELSE ! everywhere else
	341	v_ice(ji,jj) = 0._wp
	342	ENDIF
	343	!
	344	END DO
[10425]	345	CALL lbc_bdy_lnk( 'bdyice', v_ice(:,:), 'V', -1., jbdy )
[8586]	346	!
	347	END SELECT
	348	!
	349	CASE DEFAULT
	350	CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' )
	351	END SELECT
	352	!
	353	END DO
	354	!
[9905]	355	IF( ln_timing ) CALL timing_stop('bdy_ice_dyn')
	356	!
[8586]	357	END SUBROUTINE bdy_ice_dyn
	358
	359	#else
	360	!!---------------------------------------------------------------------------------
	361	!! Default option Empty module
	362	!!---------------------------------------------------------------------------------
	363	CONTAINS
	364	SUBROUTINE bdy_ice( kt ) ! Empty routine
[9927]	365	IMPLICIT NONE
	366	INTEGER, INTENT( in ) :: kt
[8586]	367	WRITE(,) 'bdy_ice: You should not have seen this print! error?', kt
	368	END SUBROUTINE bdy_ice
	369	#endif
	370
	371	!!=================================================================================
	372	END MODULE bdyice

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/OCE/BDY/bdyice.F90 @ 11083

Download in other formats: