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.

limsbc.F90 in branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 @ 4148

Last change on this file since 4148 was 4148, checked in by cetlod, 10 years ago
merge in trunk changes between r3853 and r3940 and commit the changes, see ticket #1169
Property svn:keywords set to `Id`
File size: 25.1 KB

Line
1	MODULE limsbc
2	!!======================================================================
3	!! * MODULE limsbc *
4	!! computation of the flux at the sea ice/ocean interface
5	!!======================================================================
6	!! History : - ! 2006-07 (M. Vancoppelle) LIM3 original code
7	!! 3.0 ! 2008-03 (C. Tallandier) surface module
8	!! - ! 2008-04 (C. Tallandier) split in 2 + new ice-ocean coupling
9	!! 3.3 ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
10	!! ! + simplification of the ice-ocean stress calculation
11	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
12	!! 3.5 ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
13	!!----------------------------------------------------------------------
14	#if defined key_lim3
15	!!----------------------------------------------------------------------
16	!! 'key_lim3' LIM 3.0 sea-ice model
17	!!----------------------------------------------------------------------
18	!! lim_sbc_alloc : allocate the limsbc arrays
19	!! lim_sbc_init : initialisation
20	!! lim_sbc_flx : updates mass, heat and salt fluxes at the ocean surface
21	!! lim_sbc_tau : update i- and j-stresses, and its modulus at the ocean surface
22	!!----------------------------------------------------------------------
23	USE par_oce ! ocean parameters
24	USE par_ice ! ice parameters
25	USE dom_oce ! ocean domain
26	USE sbc_ice ! Surface boundary condition: sea-ice fields
27	USE sbc_oce ! Surface boundary condition: ocean fields
28	USE phycst ! physical constants
29	USE albedo ! albedo parameters
30	USE ice ! LIM sea-ice variables
31	USE lbclnk ! ocean lateral boundary condition
32	USE in_out_manager ! I/O manager
33	USE lib_mpp ! MPP library
34	USE wrk_nemo ! work arrays
35	USE prtctl ! Print control
36	USE cpl_oasis3, ONLY : lk_cpl
37	USE oce, ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
38	USE dom_ice, ONLY : tms
39	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC lim_sbc_init ! called by ice_init
45	PUBLIC lim_sbc_flx ! called by sbc_ice_lim
46	PUBLIC lim_sbc_tau ! called by sbc_ice_lim
47
48	REAL(wp) :: epsi16 = 1.e-16_wp ! constant values
49	REAL(wp) :: rzero = 0._wp
50	REAL(wp) :: rone = 1._wp
51
52	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2]
53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmod_io ! modulus of the ice-ocean velocity [m/s]
54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soce_0 , sice_0 ! cst SSS and ice salinity (levitating sea-ice)
55
56	!! * Substitutions
57	# include "vectopt_loop_substitute.h90"
58	!!----------------------------------------------------------------------
59	!! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
60	!! $Id$
61	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
62	!!----------------------------------------------------------------------
63	CONTAINS
64
65	INTEGER FUNCTION lim_sbc_alloc()
66	!!-------------------------------------------------------------------
67	!! * ROUTINE lim_sbc_alloc *
68	!!-------------------------------------------------------------------
69	ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) , &
70	& sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=lim_sbc_alloc)
71	!
72	IF( lk_mpp ) CALL mpp_sum( lim_sbc_alloc )
73	IF( lim_sbc_alloc /= 0 ) CALL ctl_warn('lim_sbc_alloc: failed to allocate arrays')
74	END FUNCTION lim_sbc_alloc
75
76
77	SUBROUTINE lim_sbc_flx( kt )
78	!!-------------------------------------------------------------------
79	!! * ROUTINE lim_sbc_flx *
80	!!
81	!! ** Purpose : Update the surface ocean boundary condition for heat
82	!! salt and mass over areas where sea-ice is non-zero
83	!!
84	!! ** Action : - computes the heat and freshwater/salt fluxes
85	!! at the ice-ocean interface.
86	!! - Update the ocean sbc
87	!!
88	!! ** Outputs : - qsr : sea heat flux: solar
89	!! - qns : sea heat flux: non solar
90	!! - emp : freshwater budget: volume flux
91	!! - sfx : salt flux
92	!! - fr_i : ice fraction
93	!! - tn_ice : sea-ice surface temperature
94	!! - alb_ice : sea-ice alberdo (lk_cpl=T)
95	!!
96	!! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
97	!! Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
98	!!---------------------------------------------------------------------
99	INTEGER, INTENT(in) :: kt ! number of iteration
100	!
101	INTEGER :: ji, jj ! dummy loop indices
102	INTEGER :: ierr, ifvt, i1mfr, idfr ! local integer
103	INTEGER :: iflt, ial , iadv , ifral, ifrdv ! - -
104	REAL(wp) :: zinda, zemp, zemp_snow, zfmm ! local scalars
105	REAL(wp) :: zemp_snw ! - -
106	REAL(wp) :: zfcm1 , zfcm2 ! - -
107	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb, zalbp ! 2D/3D workspace
108	!!---------------------------------------------------------------------
109
110	IF( lk_cpl ) CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
111
112	!------------------------------------------!
113	! heat flux at the ocean surface !
114	!------------------------------------------!
115	! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
116	! changed to old_frld and old ht_i
117
118	DO jj = 1, jpj
119	DO ji = 1, jpi
120	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
121	ifvt = zinda * MAX( rzero , SIGN( rone, -phicif (ji,jj) ) ) !subscripts are bad here
122	i1mfr = 1.0 - MAX( rzero , SIGN( rone , - ( at_i(ji,jj) ) ) )
123	idfr = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
124	iflt = zinda * (1 - i1mfr) * (1 - ifvt )
125	ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
126	iadv = ( 1 - i1mfr ) * zinda
127	ifral = ( 1 - i1mfr * ( 1 - ial ) )
128	ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
129
130	! switch --- 1.0 ---------------- 0.0 --------------------
131	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132	! zinda \| if pfrld = 1 \| if pfrld < 1 \|
133	! -> ifvt\| if pfrld old_ht_i
134	! i1mfr \| if frld = 1 \| if frld < 1 \|
135	! idfr \| if frld <= pfrld \| if frld > pfrld \|
136	! iflt \|
137	! ial \|
138	! iadv \|
139	! ifral
140	! ifrdv
141
142	! computation the solar flux at ocean surface
143	zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj)
144	! fstric Solar flux transmitted trough the ice
145	! qsr Net short wave heat flux on free ocean
146	! new line
147	fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
148
149	! computation the non solar heat flux at ocean surface
150	zfcm2 = - zfcm1 & ! ???
151	& + iflt * fscmbq(ji,jj) & ! total ablation: heat given to the ocean
152	& + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice &
153	& + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice &
154	& + fhmec(ji,jj) & ! snow melt when ridging
155	& + fheat_mec(ji,jj) & ! ridge formation
156	& + fheat_res(ji,jj) ! residual heat flux
157	! qcmif Energy needed to bring the ocean surface layer until its freezing (ok)
158	! qldif heat balance of the lead (or of the open ocean)
159	! qfvbq latent heat uptake/release after accretion/ablation
160	! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
161
162	IF( num_sal == 2 ) zfcm2 = zfcm2 + fhbri(ji,jj) ! add contribution due to brine drainage
163
164	! bottom radiative component is sent to the computation of the oceanic heat flux
165	fsbbq(ji,jj) = ( 1._wp - ( ifvt + iflt ) ) * fscmbq(ji,jj)
166
167	! used to compute the oceanic heat flux at the next time step
168	qsr(ji,jj) = zfcm1 ! solar heat flux
169	qns(ji,jj) = zfcm2 - fdtcn(ji,jj) ! non solar heat flux
170	! ! fdtcn : turbulent oceanic heat flux
171
172	!!gm this IF prevents the vertorisation of the whole loop
173	IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN
174	WRITE(numout,*) ' lim_sbc : heat fluxes '
175	WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx)
176	WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx)
177	WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx)
178	WRITE(numout,*)
179	WRITE(numout,*) ' qns : ', qns(jiindx,jjindx)
180	WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)
181	WRITE(numout,*) ' ifral : ', ifral
182	WRITE(numout,*) ' ial : ', ial
183	WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx)
184	WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx)
185	WRITE(numout,) ' qcmif / dt: ', qcmif(jiindx,jjindx) r1_rdtice
186	WRITE(numout,) ' qldif / dt: ', qldif(jiindx,jjindx) r1_rdtice
187	WRITE(numout,*) ' ifrdv : ', ifrdv
188	WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx)
189	WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx)
190	WRITE(numout,) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) r1_rdtice
191	WRITE(numout,) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) r1_rdtice
192	WRITE(numout,*) ' '
193	WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)
194	WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx)
195	WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)
196	WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx)
197	WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
198	ENDIF
199	!!gm end
200	END DO
201	END DO
202
203	!------------------------------------------!
204	! mass flux at the ocean surface !
205	!------------------------------------------!
206
207	!!gm optimisation: this loop have to be merged with the previous one
208	DO jj = 1, jpj
209	DO ji = 1, jpi
210	! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
211	! -------------------------------------------------------------------------------------
212	! The idea of this approach is that the system that we consider is the ICE-OCEAN system
213	! Thus FW flux = External ( E-P+snow melt)
214	! Salt flux = Exchanges in the ice-ocean system then converted into FW
215	! Associated to Ice formation AND Ice melting
216	! Even if i see Ice melting as a FW and SALT flux
217	!
218
219	! computing freshwater exchanges at the ice/ocean interface
220	zemp = emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction
221	& - tprecip(ji,jj) * at_i(ji,jj) & ! all precipitation reach the ocean
222	& + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! except solid precip intercepted by sea-ice
223	& - fmmec(ji,jj) ! snow falling when ridging
224
225	! mass flux at the ocean/ice interface (sea ice fraction)
226	zemp_snw = rdm_snw(ji,jj) * r1_rdtice ! snow melting = pure water that enters the ocean
227	zfmm = rdm_ice(ji,jj) * r1_rdtice ! Freezing minus mesting
228
229	fmmflx(ji,jj) = zfmm ! F/M mass flux save at least for biogeochemical model
230
231	emp(ji,jj) = zemp + zemp_snw + zfmm ! mass flux + F/M mass flux (always ice/ocean mass exchange)
232
233	! correcting brine salt fluxes (zinda = 1 if pfrld=1 , =0 otherwise)
234	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
235	sfx_bri(ji,jj) = zinda * sfx_bri(ji,jj)
236	END DO
237	END DO
238
239	!------------------------------------------!
240	! salt flux at the ocean surface !
241	!------------------------------------------!
242
243	IF( num_sal == 2 ) THEN ! variable ice salinity: brine drainage included in the salt flux
244	sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) + sfx_bri(:,:)
245	ELSE ! constant ice salinity:
246	sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:)
247	ENDIF
248	!-----------------------------------------------!
249	! mass of snow and ice per unit area !
250	!-----------------------------------------------!
251	IF( nn_ice_embd /= 0 ) THEN ! embedded sea-ice (mass required)
252	snwice_mass_b(:,:) = snwice_mass(:,:) ! save mass from the previous ice time step
253	! ! new mass per unit area
254	snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
255	! ! time evolution of snow+ice mass
256	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
257	ENDIF
258
259	!-----------------------------------------------!
260	! Storing the transmitted variables !
261	!-----------------------------------------------!
262	fr_i (:,:) = at_i(:,:) ! Sea-ice fraction
263	tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature
264
265	!------------------------------------------------!
266	! Computation of snow/ice and ocean albedo !
267	!------------------------------------------------!
268	IF( lk_cpl ) THEN ! coupled case
269	CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb ) ! snow/ice albedo
270	!
271	alb_ice(:,:,:) = 0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:) ! Ice albedo (mean clear and overcast skys)
272	ENDIF
273
274	IF(ln_ctl) THEN
275	CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' )
276	CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=sfx , clinfo2=' sfx : ' )
277	CALL prt_ctl( tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ' )
278	CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
279	ENDIF
280	!
281	IF( lk_cpl ) CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
282	!
283	END SUBROUTINE lim_sbc_flx
284
285
286	SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
287	!!-------------------------------------------------------------------
288	!! * ROUTINE lim_sbc_tau *
289	!!
290	!! ** Purpose : Update the ocean surface stresses due to the ice
291	!!
292	!! ** Action : * at each ice time step (every nn_fsbc time step):
293	!! - compute the modulus of ice-ocean relative velocity
294	!! (rhoCd) at T-point (C-grid) or I-point (B-grid)
295	!! tmod_io = rhoco * \| U_ice-U_oce \|
296	!! - update the modulus of stress at ocean surface
297	!! taum = frld * taum + (1-frld) * tmod_io * \| U_ice-U_oce \|
298	!! * at each ocean time step (every kt):
299	!! compute linearized ice-ocean stresses as
300	!! Utau = tmod_io * \| U_ice - pU_oce \|
301	!! using instantaneous current ocean velocity (usually before)
302	!!
303	!! NB: - ice-ocean rotation angle no more allowed
304	!! - here we make an approximation: taum is only computed every ice time step
305	!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
306	!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
307	!!
308	!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
309	!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
310	!!---------------------------------------------------------------------
311	INTEGER , INTENT(in) :: kt ! ocean time-step index
312	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
313	!!
314	INTEGER :: ji, jj ! dummy loop indices
315	REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar
316	REAL(wp) :: zat_v, zvtau_ice, zv_t ! - -
317	!!---------------------------------------------------------------------
318	!
319	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
320	!CDIR NOVERRCHK
321	DO jj = 2, jpjm1 !* update the modulus of stress at ocean surface (T-point)
322	!CDIR NOVERRCHK
323	DO ji = fs_2, fs_jpim1
324	! ! 2*(U_ice-U_oce) at T-point
325	zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
326	zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
327	! ! \|U_ice-U_oce\|^2
328	zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
329	! ! update the ocean stress modulus
330	taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
331	tmod_io(ji,jj) = rhoco * SQRT( zmodt ) ! rhoco * \|U_ice-U_oce\| at T-point
332	END DO
333	END DO
334	CALL lbc_lnk( taum, 'T', 1. ) ; CALL lbc_lnk( tmod_io, 'T', 1. )
335	!
336	utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step
337	vtau_oce(:,:) = vtau(:,:)
338	!
339	ENDIF
340	!
341	! !== every ocean time-step ==!
342	!
343	DO jj = 2, jpjm1 !* update the stress WITHOUT a ice-ocean rotation angle
344	DO ji = fs_2, fs_jpim1 ! Vect. Opt.
345	zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp ! ice area at u and V-points
346	zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
347	! ! linearized quadratic drag formulation
348	zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
349	zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
350	! ! stresses at the ocean surface
351	utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
352	vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
353	END DO
354	END DO
355	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
356	!
357	IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, &
358	& tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask )
359	!
360	END SUBROUTINE lim_sbc_tau
361
362
363	SUBROUTINE lim_sbc_init
364	!!-------------------------------------------------------------------
365	!! * ROUTINE lim_sbc_init *
366	!!
367	!! ** Purpose : Preparation of the file ice_evolu for the output of
368	!! the temporal evolution of key variables
369	!!
370	!! ** input : Namelist namicedia
371	!!-------------------------------------------------------------------
372	!
373	INTEGER :: ji, jj ! dummy loop indices
374	REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar
375	IF(lwp) WRITE(numout,*)
376	IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition'
377	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ '
378
379	! ! allocate lim_sbc array
380	IF( lim_sbc_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
381	!
382	soce_0(:,:) = soce ! constant SSS and ice salinity used in levitating sea-ice case
383	sice_0(:,:) = sice
384	!
385	IF( cp_cfg == "orca" ) THEN ! decrease ocean & ice reference salinities in the Baltic sea
386	WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. &
387	& 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
388	soce_0(:,:) = 4._wp
389	sice_0(:,:) = 2._wp
390	END WHERE
391	ENDIF
392	! ! embedded sea ice
393	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
394	snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
395	snwice_mass_b(:,:) = snwice_mass(:,:)
396	ELSE
397	snwice_mass (:,:) = 0.0_wp ! no mass exchanges
398	snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges
399	ENDIF
400	IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart
401	& .NOT. ln_rstart ) THEN ! deplete the initial ssh below sea-ice area
402	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
403	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
404	!
405	! Note: Changed the initial values of sshb and sshn=> need to recompute ssh[u,v,f]_[b,n]
406	! which were previously set in domvvl
407	IF ( lk_vvl ) THEN ! Is this necessary? embd 2 should be restricted to vvl only???
408	DO jj = 1, jpjm1
409	DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible
410	zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
411	zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
412	zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
413	sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
414	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
415	sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
416	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
417	sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &
418	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
419	sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &
420	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
421	END DO
422	END DO
423	CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )
424	CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )
425	DO jj = 1, jpjm1
426	DO ji = 1, jpim1 ! NO Vector Opt.
427	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
428	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
429	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
430	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
431	END DO
432	END DO
433	CALL lbc_lnk( sshf_n, 'F', 1. )
434	ENDIF
435	ENDIF
436	!
437	END SUBROUTINE lim_sbc_init
438
439	#else
440	!!----------------------------------------------------------------------
441	!! Default option : Dummy module NO LIM 3.0 sea-ice model
442	!!----------------------------------------------------------------------
443	CONTAINS
444	SUBROUTINE lim_sbc ! Dummy routine
445	END SUBROUTINE lim_sbc
446	#endif
447
448	!!======================================================================
449	END MODULE limsbc

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

Download in other formats: