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/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 @ 3625

Last change on this file since 3625 was 3625, checked in by acc, 11 years ago
Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF
Property svn:keywords set to `Id`
File size: 24.9 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	emp(ji,jj) = zemp + zemp_snw + zfmm ! mass flux + F/M mass flux (always ice/ocean mass exchange)
230
231	! correcting brine salt fluxes (zinda = 1 if pfrld=1 , =0 otherwise)
232	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
233	sfx_bri(ji,jj) = zinda * sfx_bri(ji,jj)
234	END DO
235	END DO
236
237	!------------------------------------------!
238	! salt flux at the ocean surface !
239	!------------------------------------------!
240
241	IF( num_sal == 2 ) THEN ! variable ice salinity: brine drainage included in the salt flux
242	sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) + sfx_bri(:,:)
243	ELSE ! constant ice salinity:
244	sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:)
245	ENDIF
246	!-----------------------------------------------!
247	! mass of snow and ice per unit area !
248	!-----------------------------------------------!
249	IF( nn_ice_embd /= 0 ) THEN ! embedded sea-ice (mass required)
250	snwice_mass_b(:,:) = snwice_mass(:,:) ! save mass from the previous ice time step
251	! ! new mass per unit area
252	snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
253	! ! time evolution of snow+ice mass
254	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
255	ENDIF
256
257	!-----------------------------------------------!
258	! Storing the transmitted variables !
259	!-----------------------------------------------!
260	fr_i (:,:) = at_i(:,:) ! Sea-ice fraction
261	tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature
262
263	!------------------------------------------------!
264	! Computation of snow/ice and ocean albedo !
265	!------------------------------------------------!
266	IF( lk_cpl ) THEN ! coupled case
267	CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb ) ! snow/ice albedo
268	!
269	alb_ice(:,:,:) = 0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:) ! Ice albedo (mean clear and overcast skys)
270	ENDIF
271
272	IF(ln_ctl) THEN
273	CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' )
274	CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=sfx , clinfo2=' sfx : ' )
275	CALL prt_ctl( tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ' )
276	CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
277	ENDIF
278	!
279	IF( lk_cpl ) CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
280	!
281	END SUBROUTINE lim_sbc_flx
282
283
284	SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
285	!!-------------------------------------------------------------------
286	!! * ROUTINE lim_sbc_tau *
287	!!
288	!! ** Purpose : Update the ocean surface stresses due to the ice
289	!!
290	!! ** Action : * at each ice time step (every nn_fsbc time step):
291	!! - compute the modulus of ice-ocean relative velocity
292	!! (rhoCd) at T-point (C-grid) or I-point (B-grid)
293	!! tmod_io = rhoco * \| U_ice-U_oce \|
294	!! - update the modulus of stress at ocean surface
295	!! taum = frld * taum + (1-frld) * tmod_io * \| U_ice-U_oce \|
296	!! * at each ocean time step (every kt):
297	!! compute linearized ice-ocean stresses as
298	!! Utau = tmod_io * \| U_ice - pU_oce \|
299	!! using instantaneous current ocean velocity (usually before)
300	!!
301	!! NB: - ice-ocean rotation angle no more allowed
302	!! - here we make an approximation: taum is only computed every ice time step
303	!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
304	!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
305	!!
306	!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
307	!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
308	!!---------------------------------------------------------------------
309	INTEGER , INTENT(in) :: kt ! ocean time-step index
310	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
311	!!
312	INTEGER :: ji, jj ! dummy loop indices
313	REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar
314	REAL(wp) :: zat_v, zvtau_ice, zv_t ! - -
315	!!---------------------------------------------------------------------
316	!
317	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
318	!CDIR NOVERRCHK
319	DO jj = 2, jpjm1 !* update the modulus of stress at ocean surface (T-point)
320	!CDIR NOVERRCHK
321	DO ji = fs_2, fs_jpim1
322	! ! 2*(U_ice-U_oce) at T-point
323	zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
324	zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
325	! ! \|U_ice-U_oce\|^2
326	zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
327	! ! update the ocean stress modulus
328	taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
329	tmod_io(ji,jj) = rhoco * SQRT( zmodt ) ! rhoco * \|U_ice-U_oce\| at T-point
330	END DO
331	END DO
332	CALL lbc_lnk( taum, 'T', 1. ) ; CALL lbc_lnk( tmod_io, 'T', 1. )
333	!
334	utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step
335	vtau_oce(:,:) = vtau(:,:)
336	!
337	ENDIF
338	!
339	! !== every ocean time-step ==!
340	!
341	DO jj = 2, jpjm1 !* update the stress WITHOUT a ice-ocean rotation angle
342	DO ji = fs_2, fs_jpim1 ! Vect. Opt.
343	zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp ! ice area at u and V-points
344	zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
345	! ! linearized quadratic drag formulation
346	zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
347	zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
348	! ! stresses at the ocean surface
349	utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
350	vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
351	END DO
352	END DO
353	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
354	!
355	IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, &
356	& tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask )
357	!
358	END SUBROUTINE lim_sbc_tau
359
360
361	SUBROUTINE lim_sbc_init
362	!!-------------------------------------------------------------------
363	!! * ROUTINE lim_sbc_init *
364	!!
365	!! ** Purpose : Preparation of the file ice_evolu for the output of
366	!! the temporal evolution of key variables
367	!!
368	!! ** input : Namelist namicedia
369	!!-------------------------------------------------------------------
370	!
371	INTEGER :: ji, jj ! dummy loop indices
372	REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar
373	IF(lwp) WRITE(numout,*)
374	IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition'
375	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ '
376
377	! ! allocate lim_sbc array
378	IF( lim_sbc_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
379	!
380	soce_0(:,:) = soce ! constant SSS and ice salinity used in levitating sea-ice case
381	sice_0(:,:) = sice
382	!
383	IF( cp_cfg == "orca" ) THEN ! decrease ocean & ice reference salinities in the Baltic sea
384	WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. &
385	& 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
386	soce_0(:,:) = 4._wp
387	sice_0(:,:) = 2._wp
388	END WHERE
389	ENDIF
390	! ! embedded sea ice
391	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
392	snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
393	snwice_mass_b(:,:) = snwice_mass(:,:)
394	ELSE
395	snwice_mass (:,:) = 0.0_wp ! no mass exchanges
396	snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges
397	ENDIF
398	IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart
399	& .NOT. ln_rstart ) THEN ! deplete the initial ssh below sea-ice area
400	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
401	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
402	!
403	! Note: Changed the initial values of sshb and sshn=> need to recompute ssh[u,v,f]_[b,n]
404	! which were previously set in domvvl
405	IF ( lk_vvl ) THEN ! Is this necessary? embd 2 should be restricted to vvl only???
406	DO jj = 1, jpjm1
407	DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible
408	zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
409	zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
410	zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
411	sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
412	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
413	sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
414	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
415	sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &
416	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
417	sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &
418	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
419	END DO
420	END DO
421	CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )
422	CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )
423	DO jj = 1, jpjm1
424	DO ji = 1, jpim1 ! NO Vector Opt.
425	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
426	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
427	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
428	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
429	END DO
430	END DO
431	CALL lbc_lnk( sshf_n, 'F', 1. )
432	ENDIF
433	ENDIF
434	!
435	END SUBROUTINE lim_sbc_init
436
437	#else
438	!!----------------------------------------------------------------------
439	!! Default option : Dummy module NO LIM 3.0 sea-ice model
440	!!----------------------------------------------------------------------
441	CONTAINS
442	SUBROUTINE lim_sbc ! Dummy routine
443	END SUBROUTINE lim_sbc
444	#endif
445
446	!!======================================================================
447	END MODULE limsbc

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

Download in other formats: