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/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyice.F90 @ 11514

Last change on this file since 11514 was 11514, checked in by clem, 5 years ago
debug ponds with bdy (fraction instead of concentration)
Property svn:keywords set to `Id`
File size: 25.0 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	!
[11191]	57	INTEGER :: jbdy, ir ! BDY set index, rim index
	58	INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1)
	59	LOGICAL :: llrim0 ! indicate if rim 0 is treated
	60	LOGICAL, DIMENSION(4) :: llsend1, llrecv1 ! indicate how communications are to be carried out
[8586]	61	!!----------------------------------------------------------------------
[11042]	62	! controls
	63	IF( ln_timing ) CALL timing_start('bdy_ice_thd') ! timing
	64	IF( ln_icediachk ) CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[11501]	65	IF( ln_icediachk ) CALL ice_cons2D (0,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
[8586]	66	!
	67	CALL ice_var_glo2eqv
	68	!
[11210]	69	llsend1(:) = .false. ; llrecv1(:) = .false.
[11191]	70	DO ir = 1, 0, -1 ! treat rim 1 before rim 0
	71	IF( ir == 0 ) THEN ; llrim0 = .TRUE.
	72	ELSE ; llrim0 = .FALSE.
	73	END IF
	74	DO jbdy = 1, nb_bdy
	75	!
	76	SELECT CASE( cn_ice(jbdy) )
	77	CASE('none') ; CYCLE
	78	CASE('frs' ) ; CALL bdy_ice_frs( idx_bdy(jbdy), dta_bdy(jbdy), kt, jbdy, llrim0 )
	79	CASE DEFAULT
	80	CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' )
	81	END SELECT
	82	!
	83	END DO
[8586]	84	!
[11210]	85	! Update bdy points
	86	IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
	87	IF( nn_hls == 1 ) THEN ; llsend1(:) = .false. ; llrecv1(:) = .false. ; END IF
[11191]	88	DO jbdy = 1, nb_bdy
	89	IF( cn_ice(jbdy) == 'frs' ) THEN
	90	llsend1(:) = llsend1(:) .OR. lsend_bdyint(jbdy,1,:,ir) ! possibly every direction, T points
	91	llrecv1(:) = llrecv1(:) .OR. lrecv_bdyint(jbdy,1,:,ir) ! possibly every direction, T points
	92	END IF
	93	END DO ! jbdy
	94	IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN ! if need to send/recv in at least one direction
	95	! exchange 3d arrays
[11195]	96	CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. &
	97	& , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
	98	& , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1. &
	99	& , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
	100	! exchange 4d arrays : third dimension = 1 and then third dimension = jpk
	101	CALL lbc_lnk_multi( 'bdyice', t_s , 'T', 1., e_s , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
	102	CALL lbc_lnk_multi( 'bdyice', t_i , 'T', 1., e_i , 'T', 1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
[11067]	103	END IF
[11191]	104	END DO ! ir
[11067]	105	!
[9880]	106	CALL ice_cor( kt , 0 ) ! -- In case categories are out of bounds, do a remapping
[9885]	107	! ! i.e. inputs have not the same ice thickness distribution (set by rn_himean)
	108	! ! than the regional simulation
[9124]	109	CALL ice_var_agg(1)
	110	!
[11042]	111	! controls
[11501]	112	IF( ln_icectl ) CALL ice_prt ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) ! prints
[11042]	113	IF( ln_icediachk ) CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[11501]	114	IF( ln_icediachk ) CALL ice_cons2D (1,'bdy_ice_thd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
[11042]	115	IF( ln_timing ) CALL timing_stop ('bdy_ice_thd') ! timing
[8586]	116	!
	117	END SUBROUTINE bdy_ice
	118
	119
[11191]	120	SUBROUTINE bdy_ice_frs( idx, dta, kt, jbdy, llrim0 )
[8586]	121	!!------------------------------------------------------------------------------
	122	!! * SUBROUTINE bdy_ice_frs *
	123	!!
[9890]	124	!! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields
[8586]	125	!!
	126	!! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
	127	!! dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
	128	!!------------------------------------------------------------------------------
[11191]	129	TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices
	130	TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data
	131	INTEGER, INTENT(in) :: kt ! main time-step counter
	132	INTEGER, INTENT(in) :: jbdy ! BDY set index
	133	LOGICAL, INTENT(in) :: llrim0 ! indicate if rim 0 is treated
[8586]	134	!
	135	INTEGER :: jpbound ! 0 = incoming ice
	136	! ! 1 = outgoing ice
[11191]	137	INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1)
[9890]	138	INTEGER :: i_bdy, jgrd ! dummy loop indices
	139	INTEGER :: ji, jj, jk, jl, ib, jb
[8586]	140	REAL(wp) :: zwgt, zwgt1 ! local scalar
	141	REAL(wp) :: ztmelts, zdh
[11048]	142	REAL(wp), POINTER :: flagu, flagv ! short cuts
[8586]	143	!!------------------------------------------------------------------------------
	144	!
	145	jgrd = 1 ! Everything is at T-points here
[11191]	146	IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(jgrd)
	147	ELSE ; ibeg = idx%nblenrim0(jgrd)+1 ; iend = idx%nblenrim(jgrd)
	148	END IF
[8586]	149	!
	150	DO jl = 1, jpl
[11191]	151	DO i_bdy = ibeg, iend
[9890]	152	ji = idx%nbi(i_bdy,jgrd)
	153	jj = idx%nbj(i_bdy,jgrd)
	154	zwgt = idx%nbw(i_bdy,jgrd)
	155	zwgt1 = 1.e0 - idx%nbw(i_bdy,jgrd)
[11507]	156	a_i (ji,jj, jl) = ( a_i (ji,jj, jl) * zwgt1 + dta%a_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice concentration
	157	h_i (ji,jj, jl) = ( h_i (ji,jj, jl) * zwgt1 + dta%h_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice depth
	158	h_s (ji,jj, jl) = ( h_s (ji,jj, jl) * zwgt1 + dta%h_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow depth
	159	t_i (ji,jj,:,jl) = ( t_i (ji,jj,:,jl) * zwgt1 + dta%t_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice temperature
	160	t_s (ji,jj,:,jl) = ( t_s (ji,jj,:,jl) * zwgt1 + dta%t_s(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Snow temperature
	161	t_su(ji,jj, jl) = ( t_su(ji,jj, jl) * zwgt1 + dta%tsu(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Surf temperature
	162	s_i (ji,jj, jl) = ( s_i (ji,jj, jl) * zwgt1 + dta%s_i(i_bdy,jl) * zwgt ) * tmask(ji,jj,1) ! Ice salinity
	163	sz_i(ji,jj,:,jl) = s_i (ji,jj, jl)
	164	!
[8586]	165	! -----------------
	166	! Pathological case
	167	! -----------------
	168	! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to
	169	! very large transformation from snow to ice (see icethd_dh.F90)
	170
	171	! Then, a) transfer the snow excess into the ice (different from icethd_dh)
[9935]	172	zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )
[8586]	173	! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
[9935]	174	!zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi )
[8586]	175
	176	! recompute h_i, h_s
	177	h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh )
[9935]	178	h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )
[8586]	179
	180	ENDDO
	181	ENDDO
	182
	183	DO jl = 1, jpl
[11191]	184	DO i_bdy = ibeg, iend
[9890]	185	ji = idx%nbi(i_bdy,jgrd)
	186	jj = idx%nbj(i_bdy,jgrd)
[11048]	187	flagu => idx%flagu(i_bdy,jgrd)
	188	flagv => idx%flagv(i_bdy,jgrd)
[8586]	189	! condition on ice thickness depends on the ice velocity
	190	! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
[9890]	191	jpbound = 0 ; ib = ji ; jb = jj
[8586]	192	!
[11049]	193	IF( flagu == 1. ) THEN
[11071]	194	IF( ji+1 > jpi ) CYCLE
[11049]	195	IF( u_ice(ji ,jj ) < 0. ) jpbound = 1 ; ib = ji+1
	196	END IF
	197	IF( flagu == -1. ) THEN
[11071]	198	IF( ji-1 < 1 ) CYCLE
[11049]	199	IF( u_ice(ji-1,jj ) < 0. ) jpbound = 1 ; ib = ji-1
	200	END IF
	201	IF( flagv == 1. ) THEN
[11071]	202	IF( jj+1 > jpj ) CYCLE
[11049]	203	IF( v_ice(ji ,jj ) < 0. ) jpbound = 1 ; jb = jj+1
	204	END IF
	205	IF( flagv == -1. ) THEN
	206	IF( jj-1 < 1 ) CYCLE
	207	IF( v_ice(ji ,jj-1) < 0. ) jpbound = 1 ; jb = jj-1
	208	END IF
[8586]	209	!
[9890]	210	IF( nn_ice_dta(jbdy) == 0 ) jpbound = 0 ; ib = ji ; jb = jj ! case ice boundaries = initial conditions
	211	! ! do not make state variables dependent on velocity
[8586]	212	!
[9890]	213	IF( a_i(ib,jb,jl) > 0._wp ) THEN ! there is ice at the boundary
[8586]	214	!
[11507]	215	a_i (ji,jj, jl) = a_i (ib,jb, jl)
	216	h_i (ji,jj, jl) = h_i (ib,jb, jl)
	217	h_s (ji,jj, jl) = h_s (ib,jb, jl)
	218	t_i (ji,jj,:,jl) = t_i (ib,jb,:,jl)
	219	t_s (ji,jj,:,jl) = t_s (ib,jb,:,jl)
	220	t_su(ji,jj, jl) = t_su(ib,jb, jl)
	221	s_i (ji,jj, jl) = s_i (ib,jb, jl)
	222	sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
[8586]	223	!
[11514]	224	! ice age and ponds
	225	IF( jpbound == 0 ) THEN ! velocity is inward
	226	oa_i(ji,jj,jl) = rice_age(jbdy) * a_i(ji,jj,jl)
[11510]	227	IF( ln_pnd_CST .OR. ln_pnd_H12 ) THEN
[11514]	228	a_ip_frac(ji,jj,jl) = rice_apnd(jbdy)
	229	a_ip (ji,jj,jl) = rice_apnd(jbdy) * a_i (ji,jj,jl)
	230	v_ip (ji,jj,jl) = rice_hpnd(jbdy) * a_ip(ji,jj,jl)
[11510]	231	ENDIF
[11514]	232	ELSEIF( jpbound == 1 ) THEN ! velocity is outward
	233	oa_i(ji,jj,jl) = oa_i(ib,jb,jl)
[11510]	234	IF( ln_pnd_CST .OR. ln_pnd_H12 ) THEN
[11514]	235	a_ip_frac(ji,jj,jl) = a_ip_frac(ib,jb,jl)
	236	a_ip (ji,jj,jl) = a_ip (ib,jb,jl)
	237	v_ip (ji,jj,jl) = v_ip (ib,jb,jl)
[11510]	238	ENDIF
[11514]	239	ENDIF
[9885]	240	!
[9888]	241	IF( nn_icesal == 1 ) THEN ! if constant salinity
	242	s_i (ji,jj ,jl) = rn_icesal
	243	sz_i(ji,jj,:,jl) = rn_icesal
	244	ENDIF
[8586]	245	!
[9888]	246	! global fields
	247	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl) ! volume ice
	248	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl) ! volume snw
	249	sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
[8586]	250	DO jk = 1, nlay_s
[9935]	251	e_s(ji,jj,jk,jl) = rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) ! enthalpy in J/m3
[9888]	252	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s ! enthalpy in J/m2
	253	END DO
[8586]	254	DO jk = 1, nlay_i
[9935]	255	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) ! Melting temperature in C
[9888]	256	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
	257	!
[9935]	258	e_i(ji,jj,jk,jl) = rhoi * ( rcpi * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) ) & ! enthalpy in J/m3
	259	& + rLfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) ) &
	260	& - rcp * ztmelts )
[9888]	261	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]	262	END DO
	263	!
[9888]	264	ELSE ! no ice at the boundary
[9885]	265	!
[9888]	266	a_i (ji,jj, jl) = 0._wp
	267	h_i (ji,jj, jl) = 0._wp
	268	h_s (ji,jj, jl) = 0._wp
	269	oa_i(ji,jj, jl) = 0._wp
	270	a_ip(ji,jj, jl) = 0._wp
	271	v_ip(ji,jj, jl) = 0._wp
	272	t_su(ji,jj, jl) = rt0
	273	t_s (ji,jj,:,jl) = rt0
	274	t_i (ji,jj,:,jl) = rt0
[11514]	275
	276	a_ip_frac(ji,jj,jl) = 0._wp
	277	h_ip (ji,jj,jl) = 0._wp
	278	a_ip (ji,jj,jl) = 0._wp
	279	v_ip (ji,jj,jl) = 0._wp
[9888]	280
	281	IF( nn_icesal == 1 ) THEN ! if constant salinity
	282	s_i (ji,jj ,jl) = rn_icesal
	283	sz_i(ji,jj,:,jl) = rn_icesal
	284	ELSE ! if variable salinity
	285	s_i (ji,jj,jl) = rn_simin
	286	sz_i(ji,jj,:,jl) = rn_simin
	287	ENDIF
	288	!
	289	! global fields
	290	v_i (ji,jj, jl) = 0._wp
	291	v_s (ji,jj, jl) = 0._wp
	292	sv_i(ji,jj, jl) = 0._wp
	293	e_s (ji,jj,:,jl) = 0._wp
	294	e_i (ji,jj,:,jl) = 0._wp
	295
[8586]	296	ENDIF
[9888]	297
[8586]	298	END DO
	299	!
[9888]	300	END DO ! jl
[8586]	301	!
	302	END SUBROUTINE bdy_ice_frs
	303
	304
	305	SUBROUTINE bdy_ice_dyn( cd_type )
	306	!!------------------------------------------------------------------------------
	307	!! * SUBROUTINE bdy_ice_dyn *
	308	!!
[9890]	309	!! ** Purpose : Apply dynamics boundary conditions for sea-ice.
[8586]	310	!!
[9890]	311	!! ** Method : if this adjacent grid point is not ice free, then u_ice and v_ice take its value
	312	!! if is ice free, then u_ice and v_ice are unchanged by BDY
	313	!! they keep values calculated in rheology
	314	!!
[8586]	315	!!------------------------------------------------------------------------------
	316	CHARACTER(len=1), INTENT(in) :: cd_type ! nature of velocity grid-points
	317	!
[11191]	318	INTEGER :: i_bdy, jgrd ! dummy loop indices
	319	INTEGER :: ji, jj ! local scalar
	320	INTEGER :: jbdy, ir ! BDY set index, rim index
	321	INTEGER :: ibeg, iend ! length of rim to be treated (rim 0 or rim 1)
[8586]	322	REAL(wp) :: zmsk1, zmsk2, zflag
[11071]	323	LOGICAL, DIMENSION(4) :: llsend2, llrecv2, llsend3, llrecv3 ! indicate how communications are to be carried out
[8586]	324	!!------------------------------------------------------------------------------
[9905]	325	IF( ln_timing ) CALL timing_start('bdy_ice_dyn')
[8586]	326	!
[11210]	327	llsend2(:) = .false. ; llrecv2(:) = .false.
	328	llsend3(:) = .false. ; llrecv3(:) = .false.
[11191]	329	DO ir = 1, 0, -1
	330	DO jbdy = 1, nb_bdy
[8586]	331	!
[11191]	332	SELECT CASE( cn_ice(jbdy) )
	333	!
	334	CASE('none')
	335	CYCLE
	336	!
	337	CASE('frs')
	338	!
	339	IF( nn_ice_dta(jbdy) == 0 ) CYCLE ! case ice boundaries = initial conditions
	340	! ! do not change ice velocity (it is only computed by rheology)
	341	SELECT CASE ( cd_type )
	342	!
	343	CASE ( 'U' )
	344	jgrd = 2 ! u velocity
	345	IF( ir == 0 ) THEN ; ibeg = 1 ; iend = idx_bdy(jbdy)%nblenrim0(jgrd)
	346	ELSE ; ibeg = idx_bdy(jbdy)%nblenrim0(jgrd)+1 ; iend = idx_bdy(jbdy)%nblenrim(jgrd)
[11048]	347	END IF
[11191]	348	DO i_bdy = ibeg, iend
	349	ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
	350	jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
	351	zflag = idx_bdy(jbdy)%flagu(i_bdy,jgrd)
	352	! i-1 i i \| ! i i i+1 \| ! i i i+1 \|
	353	! > ice > \| ! o > ice \| ! o > o \|
	354	! => set at u_ice(i-1) ! => set to O ! => unchanged
	355	IF( zflag == -1. .AND. ji > 1 .AND. ji < jpi ) THEN
	356	IF ( vt_i(ji ,jj) > 0. ) THEN ; u_ice(ji,jj) = u_ice(ji-1,jj)
	357	ELSEIF( vt_i(ji+1,jj) > 0. ) THEN ; u_ice(ji,jj) = 0._wp
	358	END IF
[11048]	359	END IF
[11191]	360	! \| i i+1 i+1 ! \| i i i+1 ! \| i i i+1
	361	! \| > ice > ! \| ice > o ! \| o > o
	362	! => set at u_ice(i+1) ! => set to O ! => unchanged
	363	IF( zflag == 1. .AND. ji+1 < jpi+1 ) THEN
	364	IF ( vt_i(ji+1,jj) > 0. ) THEN ; u_ice(ji,jj) = u_ice(ji+1,jj)
	365	ELSEIF( vt_i(ji ,jj) > 0. ) THEN ; u_ice(ji,jj) = 0._wp
	366	END IF
	367	END IF
	368	!
	369	IF( zflag == 0. ) u_ice(ji,jj) = 0._wp ! u_ice = 0 if north/south bdy
	370	!
	371	END DO
	372	!
	373	CASE ( 'V' )
	374	jgrd = 3 ! v velocity
	375	IF( ir == 0 ) THEN ; ibeg = 1 ; iend = idx_bdy(jbdy)%nblenrim0(jgrd)
	376	ELSE ; ibeg = idx_bdy(jbdy)%nblenrim0(jgrd)+1 ; iend = idx_bdy(jbdy)%nblenrim(jgrd)
[11048]	377	END IF
[11191]	378	DO i_bdy = ibeg, iend
	379	ji = idx_bdy(jbdy)%nbi(i_bdy,jgrd)
	380	jj = idx_bdy(jbdy)%nbj(i_bdy,jgrd)
	381	zflag = idx_bdy(jbdy)%flagv(i_bdy,jgrd)
[11507]	382	! ! ice (jj+1) ! o (jj+1)
[11191]	383	! ^ (jj ) ! ^ (jj ) ! ^ (jj )
	384	! ice (jj ) ! o (jj ) ! o (jj )
	385	! ^ (jj-1) ! !
	386	! => set to u_ice(jj-1) ! => set to 0 ! => unchanged
	387	IF( zflag == -1. .AND. jj > 1 .AND. jj < jpj ) THEN
	388	IF ( vt_i(ji,jj ) > 0. ) THEN ; v_ice(ji,jj) = v_ice(ji,jj-1)
	389	ELSEIF( vt_i(ji,jj+1) > 0. ) THEN ; v_ice(ji,jj) = 0._wp
	390	END IF
[11048]	391	END IF
[11191]	392	! ^ (jj+1) ! !
	393	! ice (jj+1) ! o (jj+1) ! o (jj+1)
	394	! ^ (jj ) ! ^ (jj ) ! ^ (jj )
	395	! ________________ ! ____ice___(jj )_ ! _____o____(jj )
	396	! => set to u_ice(jj+1) ! => set to 0 ! => unchanged
	397	IF( zflag == 1. .AND. jj < jpj ) THEN
	398	IF ( vt_i(ji,jj+1) > 0. ) THEN ; v_ice(ji,jj) = v_ice(ji,jj+1)
	399	ELSEIF( vt_i(ji,jj ) > 0. ) THEN ; v_ice(ji,jj) = 0._wp
	400	END IF
	401	END IF
	402	!
	403	IF( zflag == 0. ) v_ice(ji,jj) = 0._wp ! v_ice = 0 if west/east bdy
	404	!
	405	END DO
[8586]	406	!
[11191]	407	END SELECT
[8586]	408	!
[11191]	409	CASE DEFAULT
	410	CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' )
[8586]	411	END SELECT
	412	!
[11191]	413	END DO ! jbdy
[8586]	414	!
[11191]	415	SELECT CASE ( cd_type )
	416	CASE ( 'U' )
[11210]	417	IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
	418	IF( nn_hls == 1 ) THEN ; llsend2(:) = .false. ; llrecv2(:) = .false. ; END IF
[11191]	419	DO jbdy = 1, nb_bdy
	420	IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN
	421	llsend2(:) = llsend2(:) .OR. lsend_bdyint(jbdy,2,:,ir) ! possibly every direction, U points
	422	llsend2(1) = llsend2(1) .OR. lsend_bdyext(jbdy,2,1,ir) ! neighbour might search point towards its west bdy
	423	llrecv2(:) = llrecv2(:) .OR. lrecv_bdyint(jbdy,2,:,ir) ! possibly every direction, U points
	424	llrecv2(2) = llrecv2(2) .OR. lrecv_bdyext(jbdy,2,2,ir) ! might search point towards east bdy
	425	END IF
	426	END DO
	427	IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN ! if need to send/recv in at least one direction
[11195]	428	CALL lbc_lnk( 'bdyice', u_ice, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
[11067]	429	END IF
[11191]	430	CASE ( 'V' )
[11210]	431	IF( nn_hls > 1 .AND. ir == 1 ) CYCLE ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
	432	IF( nn_hls == 1 ) THEN ; llsend3(:) = .false. ; llrecv3(:) = .false. ; END IF
[11191]	433	DO jbdy = 1, nb_bdy
	434	IF( cn_ice(jbdy) == 'frs' .AND. nn_ice_dta(jbdy) /= 0 ) THEN
	435	llsend3(:) = llsend3(:) .OR. lsend_bdyint(jbdy,3,:,ir) ! possibly every direction, V points
	436	llsend3(3) = llsend3(3) .OR. lsend_bdyext(jbdy,3,3,ir) ! neighbour might search point towards its south bdy
	437	llrecv3(:) = llrecv3(:) .OR. lrecv_bdyint(jbdy,3,:,ir) ! possibly every direction, V points
	438	llrecv3(4) = llrecv3(4) .OR. lrecv_bdyext(jbdy,3,4,ir) ! might search point towards north bdy
	439	END IF
	440	END DO
	441	IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN ! if need to send/recv in at least one direction
[11195]	442	CALL lbc_lnk( 'bdyice', v_ice, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
[11067]	443	END IF
[11191]	444	END SELECT
	445	END DO ! ir
[11067]	446	!
[9905]	447	IF( ln_timing ) CALL timing_stop('bdy_ice_dyn')
	448	!
[8586]	449	END SUBROUTINE bdy_ice_dyn
	450
	451	#else
	452	!!---------------------------------------------------------------------------------
	453	!! Default option Empty module
	454	!!---------------------------------------------------------------------------------
	455	CONTAINS
	456	SUBROUTINE bdy_ice( kt ) ! Empty routine
[9927]	457	IMPLICIT NONE
	458	INTEGER, INTENT( in ) :: kt
[8586]	459	WRITE(,) 'bdy_ice: You should not have seen this print! error?', kt
	460	END SUBROUTINE bdy_ice
	461	#endif
	462
	463	!!=================================================================================
	464	END MODULE bdyice

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

Download in other formats: