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

Context Navigation

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 @ 3524

Last change on this file since 3524 was 3524, checked in by gm, 12 years ago
Branch: dev_r3385_NOCS04_HAMF; #665. add USE lib_fortran when SIGN is used (TOP,OPA,LIM2&3) ; salt flux names start with sfx_ in LIM3
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) :: zfcm1 , zfcm2 ! - -
106	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb, zalbp ! 2D/3D workspace
107	!!---------------------------------------------------------------------
108
109	IF( lk_cpl ) CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
110
111	!------------------------------------------!
112	! heat flux at the ocean surface !
113	!------------------------------------------!
114	! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
115	! changed to old_frld and old ht_i
116
117	DO jj = 1, jpj
118	DO ji = 1, jpi
119	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
120	ifvt = zinda * MAX( rzero , SIGN( rone, -phicif (ji,jj) ) ) !subscripts are bad here
121	i1mfr = 1.0 - MAX( rzero , SIGN( rone , - ( at_i(ji,jj) ) ) )
122	idfr = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
123	iflt = zinda * (1 - i1mfr) * (1 - ifvt )
124	ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
125	iadv = ( 1 - i1mfr ) * zinda
126	ifral = ( 1 - i1mfr * ( 1 - ial ) )
127	ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
128
129	! switch --- 1.0 ---------------- 0.0 --------------------
130	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131	! zinda \| if pfrld = 1 \| if pfrld < 1 \|
132	! -> ifvt\| if pfrld old_ht_i
133	! i1mfr \| if frld = 1 \| if frld < 1 \|
134	! idfr \| if frld <= pfrld \| if frld > pfrld \|
135	! iflt \|
136	! ial \|
137	! iadv \|
138	! ifral
139	! ifrdv
140
141	! computation the solar flux at ocean surface
142	zfcm1 = pfrld(ji,jj) * qsr(ji,jj) + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj)
143	! fstric Solar flux transmitted trough the ice
144	! qsr Net short wave heat flux on free ocean
145	! new line
146	fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
147
148	! computation the non solar heat flux at ocean surface
149	zfcm2 = - zfcm1 & ! ???
150	& + iflt * fscmbq(ji,jj) & ! total ablation: heat given to the ocean
151	& + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice &
152	& + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice &
153	& + fhmec(ji,jj) & ! snow melt when ridging
154	& + fheat_mec(ji,jj) & ! ridge formation
155	& + fheat_res(ji,jj) ! residual heat flux
156	! qcmif Energy needed to bring the ocean surface layer until its freezing (ok)
157	! qldif heat balance of the lead (or of the open ocean)
158	! qfvbq latent heat uptake/release after accretion/ablation
159	! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
160
161	IF( num_sal == 2 ) zfcm2 = zfcm2 + fhbri(ji,jj) ! add contribution due to brine drainage
162
163	! bottom radiative component is sent to the computation of the oceanic heat flux
164	fsbbq(ji,jj) = ( 1._wp - ( ifvt + iflt ) ) * fscmbq(ji,jj)
165
166	! used to compute the oceanic heat flux at the next time step
167	qsr(ji,jj) = zfcm1 ! solar heat flux
168	qns(ji,jj) = zfcm2 - fdtcn(ji,jj) ! non solar heat flux
169	! ! fdtcn : turbulent oceanic heat flux
170
171	!!gm this IF prevents the vertorisation of the whole loop
172	IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN
173	WRITE(numout,*) ' lim_sbc : heat fluxes '
174	WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx)
175	WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx)
176	WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx)
177	WRITE(numout,*)
178	WRITE(numout,*) ' qns : ', qns(jiindx,jjindx)
179	WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)
180	WRITE(numout,*) ' ifral : ', ifral
181	WRITE(numout,*) ' ial : ', ial
182	WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx)
183	WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx)
184	WRITE(numout,) ' qcmif / dt: ', qcmif(jiindx,jjindx) r1_rdtice
185	WRITE(numout,) ' qldif / dt: ', qldif(jiindx,jjindx) r1_rdtice
186	WRITE(numout,*) ' ifrdv : ', ifrdv
187	WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx)
188	WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx)
189	WRITE(numout,) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) r1_rdtice
190	WRITE(numout,) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) r1_rdtice
191	WRITE(numout,*) ' '
192	WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)
193	WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx)
194	WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx)
195	WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx)
196	WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
197	ENDIF
198	!!gm end
199	END DO
200	END DO
201
202	!------------------------------------------!
203	! mass flux at the ocean surface !
204	!------------------------------------------!
205
206	!!gm optimisation: this loop have to be merged with the previous one
207	DO jj = 1, jpj
208	DO ji = 1, jpi
209	! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
210	! -------------------------------------------------------------------------------------
211	! The idea of this approach is that the system that we consider is the ICE-OCEAN system
212	! Thus FW flux = External ( E-P+snow melt)
213	! Salt flux = Exchanges in the ice-ocean system then converted into FW
214	! Associated to Ice formation AND Ice melting
215	! Even if i see Ice melting as a FW and SALT flux
216	!
217
218	! computing freshwater exchanges at the ice/ocean interface
219	zemp = emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction
220	& - tprecip(ji,jj) * at_i(ji,jj) & ! all precipitation reach the ocean
221	& + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! except solid precip intercepted by sea-ice
222	& - fmmec(ji,jj) ! snow falling when ridging
223
224	! mass flux at the ocean/ice interface (sea ice fraction)
225	zemp_snw = rdm_snw(ji,jj) * r1_rdtice ! snow melting = pure water that enters the ocean
226	zfmm = rdm_ice(ji,jj) * r1_rdtice ! Freezing minus mesting
227
228	emp(ji,jj) = zemp + zemp_snw + zfmm ! mass flux + F/M mass flux (always ice/ocean mass exchange)
229
230	! correcting brine salt fluxes (zinda = 1 if pfrld=1 , =0 otherwise)
231	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
232	sfx_bri(ji,jj) = zinda * sfx_bri(ji,jj)
233	END DO
234	END DO
235
236	!------------------------------------------!
237	! salt flux at the ocean surface !
238	!------------------------------------------!
239
240	IF( num_sal == 2 ) THEN ! variable ice salinity: brine drainage included in the salt flux
241	sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) + sfx_bri(:,:)
242	ELSE ! constant ice salinity:
243	sfx(:,:) = sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:)
244	ENDIF
245	!-----------------------------------------------!
246	! mass of snow and ice per unit area !
247	!-----------------------------------------------!
248	IF( nn_ice_embd /= 0 ) THEN ! embedded sea-ice (mass required)
249	snwice_mass_b(:,:) = snwice_mass(:,:) ! save mass from the previous ice time step
250	! ! new mass per unit area
251	snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
252	! ! time evolution of snow+ice mass
253	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
254	ENDIF
255
256	!-----------------------------------------------!
257	! Storing the transmitted variables !
258	!-----------------------------------------------!
259	fr_i (:,:) = at_i(:,:) ! Sea-ice fraction
260	tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature
261
262	!------------------------------------------------!
263	! Computation of snow/ice and ocean albedo !
264	!------------------------------------------------!
265	IF( lk_cpl ) THEN ! coupled case
266	CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb ) ! snow/ice albedo
267	!
268	alb_ice(:,:,:) = 0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:) ! Ice albedo (mean clear and overcast skys)
269	ENDIF
270
271	IF(ln_ctl) THEN
272	CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' )
273	CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=sfx , clinfo2=' sfx : ' )
274	CALL prt_ctl( tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ' )
275	CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
276	ENDIF
277	!
278	IF( lk_cpl ) CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
279	!
280	END SUBROUTINE lim_sbc_flx
281
282
283	SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
284	!!-------------------------------------------------------------------
285	!! * ROUTINE lim_sbc_tau *
286	!!
287	!! ** Purpose : Update the ocean surface stresses due to the ice
288	!!
289	!! ** Action : * at each ice time step (every nn_fsbc time step):
290	!! - compute the modulus of ice-ocean relative velocity
291	!! (rhoCd) at T-point (C-grid) or I-point (B-grid)
292	!! tmod_io = rhoco * \| U_ice-U_oce \|
293	!! - update the modulus of stress at ocean surface
294	!! taum = frld * taum + (1-frld) * tmod_io * \| U_ice-U_oce \|
295	!! * at each ocean time step (every kt):
296	!! compute linearized ice-ocean stresses as
297	!! Utau = tmod_io * \| U_ice - pU_oce \|
298	!! using instantaneous current ocean velocity (usually before)
299	!!
300	!! NB: - ice-ocean rotation angle no more allowed
301	!! - here we make an approximation: taum is only computed every ice time step
302	!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
303	!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
304	!!
305	!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
306	!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
307	!!---------------------------------------------------------------------
308	INTEGER , INTENT(in) :: kt ! ocean time-step index
309	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
310	!!
311	INTEGER :: ji, jj ! dummy loop indices
312	REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar
313	REAL(wp) :: zat_v, zvtau_ice, zv_t ! - -
314	!!---------------------------------------------------------------------
315	!
316	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
317	!CDIR NOVERRCHK
318	DO jj = 2, jpjm1 !* update the modulus of stress at ocean surface (T-point)
319	!CDIR NOVERRCHK
320	DO ji = fs_2, fs_jpim1
321	! ! 2*(U_ice-U_oce) at T-point
322	zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
323	zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
324	! ! \|U_ice-U_oce\|^2
325	zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
326	! ! update the ocean stress modulus
327	taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
328	tmod_io(ji,jj) = rhoco * SQRT( zmodt ) ! rhoco * \|U_ice-U_oce\| at T-point
329	END DO
330	END DO
331	CALL lbc_lnk( taum, 'T', 1. ) ; CALL lbc_lnk( tmod_io, 'T', 1. )
332	!
333	utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step
334	vtau_oce(:,:) = vtau(:,:)
335	!
336	ENDIF
337	!
338	! !== every ocean time-step ==!
339	!
340	DO jj = 2, jpjm1 !* update the stress WITHOUT a ice-ocean rotation angle
341	DO ji = fs_2, fs_jpim1 ! Vect. Opt.
342	zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp ! ice area at u and V-points
343	zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
344	! ! linearized quadratic drag formulation
345	zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
346	zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
347	! ! stresses at the ocean surface
348	utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
349	vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
350	END DO
351	END DO
352	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
353	!
354	IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, &
355	& tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask )
356	!
357	END SUBROUTINE lim_sbc_tau
358
359
360	SUBROUTINE lim_sbc_init
361	!!-------------------------------------------------------------------
362	!! * ROUTINE lim_sbc_init *
363	!!
364	!! ** Purpose : Preparation of the file ice_evolu for the output of
365	!! the temporal evolution of key variables
366	!!
367	!! ** input : Namelist namicedia
368	!!-------------------------------------------------------------------
369	!
370	INTEGER :: ji, jj ! dummy loop indices
371	REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar
372	IF(lwp) WRITE(numout,*)
373	IF(lwp) WRITE(numout,*) 'lim_sbc_init : LIM-3 sea-ice - surface boundary condition'
374	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ '
375
376	! ! allocate lim_sbc array
377	IF( lim_sbc_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
378	!
379	soce_0(:,:) = soce ! constant SSS and ice salinity used in levitating sea-ice case
380	sice_0(:,:) = sice
381	!
382	IF( cp_cfg == "orca" ) THEN ! decrease ocean & ice reference salinities in the Baltic sea
383	WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. &
384	& 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
385	soce_0(:,:) = 4._wp
386	sice_0(:,:) = 2._wp
387	END WHERE
388	ENDIF
389	! ! embedded sea ice
390	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
391	snwice_mass (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
392	snwice_mass_b(:,:) = snwice_mass(:,:)
393	ELSE
394	snwice_mass (:,:) = 0.0_wp ! no mass exchanges
395	snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges
396	ENDIF
397	IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart
398	& .NOT. ln_rstart ) THEN ! deplete the initial ssh below sea-ice area
399	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
400	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
401	!
402	! Note: Changed the initial values of sshb and sshn=> need to recompute ssh[u,v,f]_[b,n]
403	! which were previously set in domvvl
404	IF ( lk_vvl ) THEN ! Is this necessary? embd 2 should be restricted to vvl only???
405	DO jj = 1, jpjm1
406	DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible
407	zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
408	zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
409	zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
410	sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
411	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
412	sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
413	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
414	sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &
415	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
416	sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &
417	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
418	END DO
419	END DO
420	CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )
421	CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )
422	DO jj = 1, jpjm1
423	DO ji = 1, jpim1 ! NO Vector Opt.
424	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
425	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
426	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
427	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
428	END DO
429	END DO
430	CALL lbc_lnk( sshf_n, 'F', 1. )
431	ENDIF
432	ENDIF
433	!
434	END SUBROUTINE lim_sbc_init
435
436	#else
437	!!----------------------------------------------------------------------
438	!! Default option : Dummy module NO LIM 3.0 sea-ice model
439	!!----------------------------------------------------------------------
440	CONTAINS
441	SUBROUTINE lim_sbc ! Dummy routine
442	END SUBROUTINE lim_sbc
443	#endif
444
445	!!======================================================================
446	END MODULE limsbc

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

Download in other formats: