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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 @ 2370

Last change on this file since 2370 was 2370, checked in by gm, 13 years ago
v3.3beta: ice-ocean stress at kt with VP & EVP (LIM-2 and -3)
Property svn:keywords set to `Id`
File size: 21.3 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	!!----------------------------------------------------------------------
12	#if defined key_lim3
13	!!----------------------------------------------------------------------
14	!! 'key_lim3' LIM 3.0 sea-ice model
15	!!----------------------------------------------------------------------
16	!! lim_sbc_flx : updates mass, heat and salt fluxes at the ocean surface
17	!! lim_sbc_tau : update i- and j-stresses, and its modulus at the ocean surface
18	!!----------------------------------------------------------------------
19	USE par_oce ! ocean parameters
20	USE par_ice ! ice parameters
21	USE dom_oce ! ocean domain
22	USE sbc_ice ! Surface boundary condition: sea-ice fields
23	USE sbc_oce ! Surface boundary condition: ocean fields
24	USE phycst ! physical constants
25	USE ice ! LIM sea-ice variables
26	USE iceini ! ???
27
28	USE lbclnk ! ocean lateral boundary condition
29	USE in_out_manager ! I/O manager
30	USE albedo ! albedo parameters
31	USE prtctl ! Print control
32
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC lim_sbc_flx ! called by sbc_ice_lim
37	PUBLIC lim_sbc_tau ! called by sbc_ice_lim
38
39	REAL(wp) :: r1_rdtice ! = 1. / rdt_ice
40	REAL(wp) :: epsi16 = 1.e-16_wp ! constant values
41	REAL(wp) :: rzero = 0._wp
42	REAL(wp) :: rone = 1._wp
43
44	REAL(wp), DIMENSION(jpi,jpj) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2]
45	REAL(wp), DIMENSION(jpi,jpj) :: tmod_io ! modulus of the ice-ocean relative velocity [m/s]
46
47	REAL(wp), DIMENSION(jpi,jpj) :: soce_0, sice_0 ! constant SSS and ice salinity used in levitating sea-ice case
48
49	!! * Substitutions
50	# include "vectopt_loop_substitute.h90"
51	!!----------------------------------------------------------------------
52	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
53	!! $Id$
54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
55	!!----------------------------------------------------------------------
56	CONTAINS
57
58	SUBROUTINE lim_sbc_flx( kt )
59	!!-------------------------------------------------------------------
60	!! * ROUTINE lim_sbc_flx *
61	!!
62	!! ** Purpose : Update the surface ocean boundary condition for heat
63	!! salt and mass over areas where sea-ice is non-zero
64	!!
65	!! ** Action : - computes the heat and freshwater/salt fluxes
66	!! at the ice-ocean interface.
67	!! - Update the ocean sbc
68	!!
69	!! ** Outputs : - qsr : sea heat flux: solar
70	!! - qns : sea heat flux: non solar
71	!! - emp : freshwater budget: volume flux
72	!! - emps : freshwater budget: concentration/dillution
73	!! - fr_i : ice fraction
74	!! - tn_ice : sea-ice surface temperature
75	!! - alb_ice : sea-ice alberdo (lk_cpl=T)
76	!!
77	!! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
78	!! Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
79	!!---------------------------------------------------------------------
80	INTEGER, INTENT(in) :: kt ! number of iteration
81	!!
82	INTEGER :: ji, jj ! dummy loop indices
83	INTEGER :: ifvt, i1mfr, idfr ! some switches
84	INTEGER :: iflt, ial, iadv, ifral, ifrdv
85	REAL(wp) :: zinda ! switch for testing the values of ice concentration
86	REAL(wp) :: zfons ! salt exchanges at the ice/ocean interface
87	REAL(wp) :: zpme ! freshwater exchanges at the ice/ocean interface
88	REAL(wp), DIMENSION(jpi,jpj) :: zfcm1 , zfcm2 ! solar/non solar heat fluxes
89	#if defined key_coupled
90	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalb ! albedo of ice under overcast sky
91	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zalbp ! albedo of ice under clear sky
92	#endif
93	!!---------------------------------------------------------------------
94
95	IF( kt == nit000 ) THEN
96	IF(lwp) WRITE(numout,*)
97	IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes'
98	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
99	!
100	r1_rdtice = 1. / rdt_ice
101	!
102	soce_0(:,:) = soce
103	sice_0(:,:) = sice
104	!
105	IF( cp_cfg == "orca" ) THEN ! decrease ocean & ice reference salinities in the Baltic sea
106	WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. &
107	& 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
108	soce_0(:,:) = 4._wp
109	sice_0(:,:) = 2._wp
110	END WHERE
111	ENDIF
112	!
113	ENDIF
114
115	!------------------------------------------!
116	! heat flux at the ocean surface !
117	!------------------------------------------!
118	! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
119	! changed to old_frld and old ht_i
120
121	DO jj = 1, jpj
122	DO ji = 1, jpi
123	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
124	ifvt = zinda * MAX( rzero , SIGN( rone, -phicif (ji,jj) ) ) !subscripts are bad here
125	i1mfr = 1.0 - MAX( rzero , SIGN( rone , - ( at_i(ji,jj) ) ) )
126	idfr = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) )
127	iflt = zinda * (1 - i1mfr) * (1 - ifvt )
128	ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
129	iadv = ( 1 - i1mfr ) * zinda
130	ifral = ( 1 - i1mfr * ( 1 - ial ) )
131	ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
132
133	! switch --- 1.0 ---------------- 0.0 --------------------
134	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135	! zinda \| if pfrld = 1 \| if pfrld < 1 \|
136	! -> ifvt\| if pfrld old_ht_i
137	! i1mfr \| if frld = 1 \| if frld < 1 \|
138	! idfr \| if frld <= pfrld \| if frld > pfrld \|
139	! iflt \|
140	! ial \|
141	! iadv \|
142	! ifral
143	! ifrdv
144
145	! computation the solar flux at ocean surface
146	zfcm1(ji,jj) = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
147	! fstric Solar flux transmitted trough the ice
148	! qsr Net short wave heat flux on free ocean
149	! new line
150	fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj)
151
152	! computation the non solar heat flux at ocean surface
153	zfcm2(ji,jj) = - zfcm1(ji,jj) &
154	& + iflt * ( fscmbq(ji,jj) ) & ! total abl -> fscmbq is given to the ocean
155	! fscmbq and ffltbif are obsolete
156	! & + iflt * ffltbif(ji,jj) !!! only if one category is used
157	& + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice &
158	& + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice &
159	& + fhmec(ji,jj) & ! new contribution due to snow melt in ridging!!
160	& + fheat_rpo(ji,jj) & ! contribution from ridge formation
161	& + fheat_res(ji,jj)
162	! fscmbq Part of the solar radiation transmitted through the ice and going to the ocean
163	! computed in limthd_zdf.F90
164	! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t
165	! qcmif Energy needed to bring the ocean surface layer until its freezing (ok)
166	! qldif heat balance of the lead (or of the open ocean)
167	! qfvbq i think this is wrong!
168	! ---> Array used to store energy in case of total lateral ablation
169	! qfvbq latent heat uptake/release after accretion/ablation
170	! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
171
172	IF ( num_sal .EQ. 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + &
173	fhbri(ji,jj) ! new contribution due to brine drainage
174
175	! bottom radiative component is sent to the computation of the
176	! oceanic heat flux
177	fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)
178
179	! used to compute the oceanic heat flux at the next time step
180	qsr(ji,jj) = zfcm1(ji,jj) ! solar heat flux
181	qns(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj) ! non solar heat flux
182	! ! fdtcn : turbulent oceanic heat flux
183
184	!!gm this IF prevents the vertorisation of the whole loop
185	IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN
186	WRITE(numout,*) ' lim_sbc : heat fluxes '
187	WRITE(numout,*) ' qsr : ', qsr(jiindx,jjindx)
188	WRITE(numout,*) ' zfcm1 : ', zfcm1(jiindx,jjindx)
189	WRITE(numout,*) ' pfrld : ', pfrld(jiindx,jjindx)
190	WRITE(numout,*) ' fstric : ', fstric (jiindx,jjindx)
191	WRITE(numout,*)
192	WRITE(numout,*) ' qns : ', qns(jiindx,jjindx)
193	WRITE(numout,*) ' zfcm2 : ', zfcm2(jiindx,jjindx)
194	WRITE(numout,*) ' zfcm1 : ', zfcm1(jiindx,jjindx)
195	WRITE(numout,*) ' ifral : ', ifral
196	WRITE(numout,*) ' ial : ', ial
197	WRITE(numout,*) ' qcmif : ', qcmif(jiindx,jjindx)
198	WRITE(numout,*) ' qldif : ', qldif(jiindx,jjindx)
199	WRITE(numout,) ' qcmif / dt: ', qcmif(jiindx,jjindx) r1_rdtice
200	WRITE(numout,) ' qldif / dt: ', qldif(jiindx,jjindx) r1_rdtice
201	WRITE(numout,*) ' ifrdv : ', ifrdv
202	WRITE(numout,*) ' qfvbq : ', qfvbq(jiindx,jjindx)
203	WRITE(numout,*) ' qdtcn : ', qdtcn(jiindx,jjindx)
204	WRITE(numout,) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) r1_rdtice
205	WRITE(numout,) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) r1_rdtice
206	WRITE(numout,*) ' '
207	WRITE(numout,*) ' fdtcn : ', fdtcn(jiindx,jjindx)
208	WRITE(numout,*) ' fhmec : ', fhmec(jiindx,jjindx)
209	WRITE(numout,*) ' fheat_rpo : ', fheat_rpo(jiindx,jjindx)
210	WRITE(numout,*) ' fhbri : ', fhbri(jiindx,jjindx)
211	WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
212	ENDIF
213	!!gm end
214	END DO
215	END DO
216
217	!------------------------------------------!
218	! mass flux at the ocean surface !
219	!------------------------------------------!
220
221	!!gm optimisation: this loop have to be merged with the previous one
222	DO jj = 1, jpj
223	DO ji = 1, jpi
224	! case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
225	! -------------------------------------------------------------------------------------
226	! The idea of this approach is that the system that we consider is the ICE-OCEAN system
227	! Thus FW flux = External ( E-P+snow melt)
228	! Salt flux = Exchanges in the ice-ocean system then converted into FW
229	! Associated to Ice formation AND Ice melting
230	! Even if i see Ice melting as a FW and SALT flux
231	!
232
233	! computing freshwater exchanges at the ice/ocean interface
234	zpme = - emp(ji,jj) * ( 1.0 - at_i(ji,jj) ) & ! evaporation over oceanic fraction
235	& + tprecip(ji,jj) * at_i(ji,jj) & ! total precipitation
236	! old fashioned way
237	! & - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) ) & ! remov. snow precip over ice
238	& - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) ) & ! remov. snow precip over ice
239	& - rdmsnif(ji,jj) * r1_rdtice & ! freshwaterflux due to snow melting
240	! new contribution from snow falling when ridging
241	& + fmmec(ji,jj)
242
243	! computing salt exchanges at the ice/ocean interface
244	! sice should be the same as computed with the ice model
245	zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice
246	! SOCE
247	zfons = ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice
248
249	!CT useless ! salt flux for constant salinity
250	!CT useless fsalt(ji,jj) = zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj)
251	! salt flux for variable salinity
252	zinda = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
253	! correcting brine and salt fluxes
254	fsbri(ji,jj) = zinda*fsbri(ji,jj)
255	! converting the salt fluxes from ice to a freshwater flux from ocean
256	fsalt_res(ji,jj) = fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )
257	fseqv(ji,jj) = fseqv(ji,jj) / ( sss_m(ji,jj) + epsi16 )
258	fsbri(ji,jj) = fsbri(ji,jj) / ( sss_m(ji,jj) + epsi16 )
259	fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )
260
261	! freshwater mass exchange (positive to the ice, negative for the ocean ?)
262	! actually it's a salt flux (so it's minus freshwater flux)
263	! if sea ice grows, zfons is positive, fsalt also
264	! POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN
265	! POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1]
266
267	emp(ji,jj) = - zpme
268	END DO
269	END DO
270
271	IF( num_sal == 2 ) THEN ! variable ice salinity: brine drainage included in the salt flux
272	emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
273	ELSE ! constant ice salinity:
274	emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
275	ENDIF
276
277	!-----------------------------------------------!
278	! Storing the transmitted variables !
279	!-----------------------------------------------!
280
281	fr_i (:,:) = at_i(:,:) ! Sea-ice fraction
282	tn_ice(:,:,:) = t_su(:,:,:) ! Ice surface temperature
283
284	#if defined key_coupled
285	!------------------------------------------------!
286	! Computation of snow/ice and ocean albedo !
287	!------------------------------------------------!
288	zalb (:,:,:) = 0.e0
289	zalbp (:,:,:) = 0.e0
290
291	CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )
292
293	alb_ice(:,:,:) = 0.5 * zalbp(:,:,:) + 0.5 * zalb (:,:,:) ! Ice albedo (mean clear and overcast skys)
294	#endif
295
296	IF(ln_ctl) THEN
297	CALL prt_ctl( tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ' )
298	CALL prt_ctl( tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=emps, clinfo2=' emps : ' )
299	CALL prt_ctl( tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ' )
300	CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
301	ENDIF
302	!
303	END SUBROUTINE lim_sbc_flx
304
305
306	SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
307	!!-------------------------------------------------------------------
308	!! * ROUTINE lim_sbc_tau *
309	!!
310	!! ** Purpose : Update the ocean surface stresses due to the ice
311	!!
312	!! ** Action : * at each ice time step (every nn_fsbc time step):
313	!! - compute the modulus of ice-ocean relative velocity
314	!! (rhoCd) at T-point (C-grid) or I-point (B-grid)
315	!! tmod_io = rhoco * \| U_ice-U_oce \|
316	!! - update the modulus of stress at ocean surface
317	!! taum = frld * taum + (1-frld) * tmod_io * \| U_ice-U_oce \|
318	!! * at each ocean time step (every kt):
319	!! compute linearized ice-ocean stresses as
320	!! Utau = tmod_io * \| U_ice - pU_oce \|
321	!! using instantaneous current ocean velocity (usually before)
322	!!
323	!! NB: - ice-ocean rotation angle no more allowed
324	!! - here we make an approximation: taum is only computed every ice time step
325	!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
326	!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
327	!!
328	!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
329	!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
330	!!---------------------------------------------------------------------
331	INTEGER , INTENT(in) :: kt ! ocean time-step index
332	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
333	!!
334	INTEGER :: ji, jj ! dummy loop indices
335	REAL(wp) :: zat_u, zutau_ice, zu_t, zmodt ! local scalar
336	REAL(wp) :: zat_v, zvtau_ice, zv_t ! - -
337	!!---------------------------------------------------------------------
338
339	IF( kt == nit000 ) THEN
340	IF(lwp) WRITE(numout,*)
341	IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM-3 sea-ice - surface ocean momentum fluxes'
342	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
343	ENDIF
344
345	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
346	!CDIR NOVERRCHK
347	DO jj = 2, jpjm1 !* update the modulus of stress at ocean surface (T-point)
348	!CDIR NOVERRCHK
349	DO ji = fs_2, fs_jpim1
350	! ! 2*(U_ice-U_oce) at T-point
351	zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
352	zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
353	! ! \|U_ice-U_oce\|^2
354	zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t )
355	! ! update the ocean stress modulus
356	taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
357	tmod_io(ji,jj) = rhoco * SQRT( zmodt ) ! rhoco * \|U_ice-U_oce\| at T-point
358	END DO
359	END DO
360	CALL lbc_lnk( taum, 'T', 1. ) ; CALL lbc_lnk( tmod_io, 'T', 1. )
361	!
362	utau_oce(:,:) = utau(:,:) !* save the air-ocean stresses at ice time-step
363	vtau_oce(:,:) = vtau(:,:)
364	!
365	ENDIF
366	!
367	! !== every ocean time-step ==!
368	!
369	DO jj = 2, jpjm1 !* update the stress WITHOUT a ice-ocean rotation angle
370	DO ji = fs_2, fs_jpim1 ! Vect. Opt.
371	zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp ! ice area at u and V-points
372	zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
373	! ! linearized quadratic drag formulation
374	zutau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
375	zvtau_ice = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
376	! ! stresses at the ocean surface
377	utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
378	vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
379	END DO
380	END DO
381	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
382	!
383	IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, &
384	& tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask )
385	!
386	END SUBROUTINE lim_sbc_tau
387
388	#else
389	!!----------------------------------------------------------------------
390	!! Default option : Dummy module NO LIM 3.0 sea-ice model
391	!!----------------------------------------------------------------------
392	CONTAINS
393	SUBROUTINE lim_sbc ! Dummy routine
394	END SUBROUTINE lim_sbc
395	#endif
396
397	!!======================================================================
398	END MODULE limsbc

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

Download in other formats: