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_2.F90 in branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago
changes after review
Property svn:keywords set to `Id`
File size: 32.7 KB

Line
1	MODULE limsbc_2
2	!!======================================================================
3	!! * MODULE limsbc_2 *
4	!! LIM-2 : updates the fluxes at the ocean surface with ice-ocean fluxes
5	!!======================================================================
6	!! History : LIM ! 2000-01 (H. Goosse) Original code
7	!! 1.0 ! 2002-07 (C. Ethe, G. Madec) re-writing F90
8	!! 3.0 ! 2006-07 (G. Madec) surface module
9	!! 3.3 ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
10	!! - ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step
11	!! 3.3.1 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation
12	!! 3.5 ! 2012-11 ((G. Madec, Y. Aksenov, A. Coward) salt and heat fluxes associated with e-p
13	!!----------------------------------------------------------------------
14	#if defined key_lim2
15	!!----------------------------------------------------------------------
16	!! 'key_lim2' LIM 2.0 sea-ice model
17	!!----------------------------------------------------------------------
18	!! lim_sbc_alloc_2 : allocate the limsbc arrays
19	!! lim_sbc_init : initialisation
20	!! lim_sbc_flx_2 : update mass, heat and salt fluxes at the ocean surface
21	!! lim_sbc_tau_2 : update i- and j-stresses, and its modulus at the ocean surface
22	!!----------------------------------------------------------------------
23	USE par_oce ! ocean parameters
24	USE phycst ! physical constants
25	USE dom_oce ! ocean domain
26	USE domvvl ! ocean vertical scale factors
27	USE dom_ice_2 ! LIM-2: ice domain
28	USE ice_2 ! LIM-2: ice variables
29	USE sbc_ice ! surface boundary condition: ice
30	USE sbc_oce ! surface boundary condition: ocean
31	USE sbccpl ! surface boundary condition: coupled interface
32	USE oce , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
33	USE albedo ! albedo parameters
34	!
35	USE lbclnk ! ocean lateral boundary condition - MPP exchanges
36	USE lib_mpp ! MPP library
37	USE wrk_nemo ! work arrays
38	USE in_out_manager ! I/O manager
39	USE iom ! I/O library
40	USE prtctl ! Print control
41	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
42
43	IMPLICIT NONE
44	PRIVATE
45
46	PUBLIC lim_sbc_init_2 ! called by ice_init_2
47	PUBLIC lim_sbc_flx_2 ! called by sbc_ice_lim_2
48	PUBLIC lim_sbc_tau_2 ! called by sbc_ice_lim_2
49
50	REAL(wp) :: r1_rdtice ! = 1. / rdt_ice
51	REAL(wp) :: epsi16 = 1.e-16_wp ! constant values
52	REAL(wp) :: rzero = 0._wp ! - -
53	REAL(wp) :: rone = 1._wp ! - -
54	!
55	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: soce_0, sice_0 ! fix SSS and ice salinity used in levitating case 0
56	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: utau_oce, vtau_oce ! air-ocean surface i- & j-stress [N/m2]
57	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmod_io ! modulus of the ice-ocean relative velocity [m/s]
58
59	!! * Substitutions
60	# include "vectopt_loop_substitute.h90"
61	!!----------------------------------------------------------------------
62	!! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2011)
63	!! $Id$
64	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
65	!!----------------------------------------------------------------------
66	CONTAINS
67
68	INTEGER FUNCTION lim_sbc_alloc_2()
69	!!-------------------------------------------------------------------
70	!! * ROUTINE lim_sbc_alloc_2 *
71	!!-------------------------------------------------------------------
72	ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) , &
73	& sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=lim_sbc_alloc_2)
74	!
75	IF( lk_mpp ) CALL mpp_sum( lim_sbc_alloc_2 )
76	IF( lim_sbc_alloc_2 /= 0 ) CALL ctl_warn('lim_sbc_alloc_2: failed to allocate arrays.')
77	!
78	END FUNCTION lim_sbc_alloc_2
79
80
81	SUBROUTINE lim_sbc_flx_2( kt )
82	!!-------------------------------------------------------------------
83	!! * ROUTINE lim_sbc_2 *
84	!!
85	!! ** Purpose : Update surface ocean boundary condition over areas
86	!! that are at least partially covered by sea-ice
87	!!
88	!! ** Action : - comput. of the momentum, heat and freshwater/salt
89	!! fluxes at the ice-ocean interface.
90	!! - Update the fluxes provided to the ocean
91	!!
92	!! ** Outputs : - qsr : sea heat flux : solar
93	!! - qns : sea heat flux : non solar (including heat content of the mass flux)
94	!! - emp : freshwater budget: mass flux
95	!! - sfx : freshwater budget: salt flux due to Freezing/Melting
96	!! - fr_i : ice fraction
97	!! - tn_ice : sea-ice surface temperature
98	!! - alb_ice : sea-ice albedo (ln_cpl=T)
99	!!
100	!! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
101	!! Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
102	!!---------------------------------------------------------------------
103	INTEGER, INTENT(in) :: kt ! number of iteration
104	!
105	INTEGER :: ji, jj ! dummy loop indices
106	INTEGER :: ii0, ii1, ij0, ij1 ! local integers
107	INTEGER :: ifvt, i1mfr, idfr, iflt ! - -
108	INTEGER :: ial, iadv, ifral, ifrdv ! - -
109	REAL(wp) :: zqsr, zqns, zfmm ! local scalars
110	REAL(wp) :: zinda, zfsalt, zemp ! - -
111	REAL(wp) :: zemp_snw, zqhc, zcd ! - -
112	REAL(wp) :: zswitch ! - -
113	REAL(wp), POINTER, DIMENSION(:,:) :: zqnsoce ! 2D workspace
114	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb, zalbp ! 2D/3D workspace
115	!!---------------------------------------------------------------------
116	!
117	CALL wrk_alloc( jpi, jpj, zqnsoce )
118	CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp )
119	!
120	SELECT CASE( nn_ice_embd ) ! levitating or embedded sea-ice option
121	CASE( 0 ) ; zswitch = 1 ! (0) old levitating sea-ice : salt exchange only
122	CASE( 1, 2 ) ; zswitch = 0 ! (1) levitating sea-ice: salt and volume exchange but no pressure effect
123	! ! (2) embedded sea-ice : salt and volume fluxes and pressure
124	END SELECT
125
126	!------------------------------------------!
127	! heat flux at the ocean surface !
128	!------------------------------------------!
129
130	!$OMP PARALLEL
131	!$OMP DO schedule(static) private(jj, ji)
132	DO jj = 1, jpj
133	DO ji = 1, jpi
134	zqnsoce(ji,jj) = qns(ji,jj)
135	END DO
136	END DO
137	!$OMP END DO NOWAIT
138
139	!$OMP DO schedule(static) private(jj,ji,zinda,ifvt,i1mfr,idfr,iflt,ial,iadv,ifral,ifrdv,zqsr,zqns,zqhc,zemp,zemp_snw,zfmm,zfsalt,zcd)
140	DO jj = 1, jpj
141	DO ji = 1, jpi
142	zinda = 1.0 - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) ) ) )
143	ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) )
144	i1mfr = 1.0 - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) )
145	idfr = 1.0 - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
146	iflt = zinda * (1 - i1mfr) * (1 - ifvt )
147	ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
148	iadv = ( 1 - i1mfr ) * zinda
149	ifral = ( 1 - i1mfr * ( 1 - ial ) )
150	ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
151
152	!!$ attempt to explain the tricky flags set above....
153	!!$ zinda = 1.0 - AINT( pfrld(ji,jj) ) ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo)
154	!!$ i1mfr = 1.0 - AINT( frld(ji,jj) ) ! = 0. if ice-free ocean else 1. (after ice thermo)
155	!!$
156	!!$ IF( phicif(ji,jj) <= 0. ) THEN ; ifvt = zinda ! = zinda if previous thermodynamic step overmelted the ice???
157	!!$ ELSE ; ifvt = 0. !
158	!!$ ENDIF
159	!!$
160	!!$ IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN ; idfr = 0. ! = 0. if lead fraction increases due to ice thermodynamics
161	!!$ ELSE ; idfr = 1.
162	!!$ ENDIF
163	!!$
164	!!$ iflt = zinda * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous time and ice-free ocean currently
165	!!$
166	!!$ ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
167	!!$ = i1mfr if ifvt = 1 i.e.
168	!!$ = idfr if ifvt = 0
169	!!$! snow no ice ice ice or nothing lead fraction increases
170	!!$! at previous now at previous
171	!!$! -> ice area increases ??? -> ice area decreases ???
172	!!$
173	!!$ iadv = ( 1 - i1mfr ) * zinda
174	!!$! pure ocean ice at
175	!!$! at current previous
176	!!$! -> = 1. if ice disapear between previous and current
177	!!$
178	!!$ ifral = ( 1 - i1mfr * ( 1 - ial ) )
179	!!$! ice at ???
180	!!$! current
181	!!$! -> ???
182	!!$
183	!!$ ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
184	!!$! ice disapear
185	!!$
186	!!$
187
188	! computation the solar flux at ocean surface
189	IF( ln_cpl ) THEN
190	zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
191	ELSE
192	zqsr = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
193	ENDIF
194	! computation the non solar heat flux at ocean surface
195	zqns = - ( 1. - thcm(ji,jj) ) * zqsr & ! part of the solar energy used in leads
196	& + iflt * ( fscmbq(ji,jj) + ffltbif(ji,jj) ) &
197	& + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice &
198	& + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice
199
200	fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj) ! store residual heat flux (to put into the ocean at the next time-step)
201	zqhc = ( rdq_snw(ji,jj) &
202	& + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice ! heat flux due to snow ( & ice heat content,
203	! ! if ice/ocean mass exchange active)
204	qsr (ji,jj) = zqsr ! solar heat flux
205	qns (ji,jj) = zqns - fdtcn(ji,jj) + zqhc ! non solar heat flux
206	!
207	! !------------------------------------------!
208	! ! mass and salt flux at the ocean surface !
209	! !------------------------------------------!
210	!
211	! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
212	! ! coupled mode:
213	IF( ln_cpl ) THEN
214	zemp = + emp_tot(ji,jj) & ! net mass flux over the grid cell (ice+ocean area)
215	& - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) ) ! minus the mass flux intercepted by sea-ice
216	ELSE
217	! ! forced mode:
218	zemp = + emp(ji,jj) * frld(ji,jj) & ! mass flux over open ocean fraction
219	& - tprecip(ji,jj) * ( 1. - frld(ji,jj) ) & ! liquid precip. over ice reaches directly the ocean
220	& + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) ) ! snow is intercepted by sea-ice (previous frld)
221	ENDIF
222	!
223	! mass flux at the ocean/ice interface (sea ice fraction)
224	zemp_snw = rdm_snw(ji,jj) * r1_rdtice ! snow melting = pure water that enters the ocean
225	zfmm = rdm_ice(ji,jj) * r1_rdtice ! Freezing minus Melting (F-M)
226
227	fmmflx(ji,jj) = zfmm ! F/M mass flux save at least for biogeochemical model
228
229	! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
230	zfsalt = - sice_0(ji,jj) * zfmm ! F-M salt exchange
231	zcd = soce_0(ji,jj) * zfmm ! concentration/dilution term due to F-M
232	!
233	! salt flux only : add concentration dilution term in salt flux and no F-M term in volume flux
234	! salt and mass fluxes : non concentration dilution term in salt flux and add F-M term in volume flux
235	sfx (ji,jj) = zfsalt + zswitch * zcd ! salt flux (+ C/D if no ice/ocean mass exchange)
236	emp (ji,jj) = zemp + zemp_snw + ( 1.- zswitch) * zfmm ! mass flux (+ F/M mass flux if ice/ocean mass exchange)
237	!
238	END DO
239	END DO
240	!$OMP END DO NOWAIT
241	!$OMP END PARALLEL
242	! !------------------------------------------!
243	! ! mass of snow and ice per unit area !
244	! !------------------------------------------!
245	IF( nn_ice_embd /= 0 ) THEN ! embedded sea-ice (mass required)
246	!$OMP PARALLEL DO schedule(static) private(jj, ji)
247	DO jj = 1, jpj
248	DO ji = 1, jpi
249	snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step
250	! ! new mass per unit area
251	snwice_mass (ji,jj) = tms(ji,jj) * ( rhosn * hsnif(ji,jj) + rhoic * hicif(ji,jj) ) * ( 1.0 - frld(ji,jj) )
252	! ! time evolution of snow+ice mass
253	snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) / rdt_ice
254	END DO
255	END DO
256	ENDIF
257
258	IF( iom_use('hflx_ice_cea' ) ) CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )
259	IF( iom_use('qns_io_cea' ) ) CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )
260	IF( iom_use('qsr_io_cea' ) ) CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
261
262	IF( iom_use('isnwmlt_cea' ) ) CALL iom_put( 'isnwmlt_cea' , rdm_snw(:,:) * r1_rdtice )
263	IF( iom_use('fsal_virt_cea') ) CALL iom_put( 'fsal_virt_cea', soce_0(:,:) * rdm_ice(:,:) * r1_rdtice )
264	IF( iom_use('fsal_real_cea') ) CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice )
265
266	!-----------------------------------------------!
267	! Coupling variables !
268	!-----------------------------------------------!
269
270	IF( ln_cpl) THEN
271	!$OMP PARALLEL DO schedule(static) private(jj, ji)
272	DO jj = 1, jpj
273	DO ji = 1, jpi
274	tn_ice(ji,jj,1) = sist(ji,jj) ! sea-ice surface temperature
275	ht_i(ji,jj,1) = hicif(ji,jj)
276	ht_s(ji,jj,1) = hsnif(ji,jj)
277	a_i(ji,jj,1) = fr_i(ji,jj)
278	END DO
279	END DO
280	! ! Computation of snow/ice and ocean albedo
281	CALL albedo_ice( tn_ice, ht_i, ht_s, zalbp, zalb )
282	!$OMP PARALLEL DO schedule(static) private(jj, ji)
283	DO jj = 1, jpj
284	DO ji = 1, jpi
285	alb_ice(ji,jj,1) = 0.5 * ( zalbp(ji,jj,1) + zalb (ji,jj,1) ) ! Ice albedo (mean clear and overcast skys)
286	END DO
287	END DO
288	IF( iom_use('icealb_cea' ) ) CALL iom_put( 'icealb_cea', alb_ice(:,:,1) * fr_i(:,:) ) ! ice albedo
289	ENDIF
290
291	IF(ln_ctl) THEN ! control print
292	CALL prt_ctl(tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ')
293	CALL prt_ctl(tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=sfx , clinfo2=' sfx : ')
294	CALL prt_ctl(tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice : ')
295	ENDIF
296	!
297	CALL wrk_dealloc( jpi, jpj, zqnsoce )
298	CALL wrk_dealloc( jpi, jpj, 1, zalb, zalbp )
299	!
300	END SUBROUTINE lim_sbc_flx_2
301
302
303	SUBROUTINE lim_sbc_tau_2( kt , pu_oce, pv_oce )
304	!!-------------------------------------------------------------------
305	!! * ROUTINE lim_sbc_tau_2 *
306	!!
307	!! ** Purpose : Update the ocean surface stresses due to the ice
308	!!
309	!! ** Action : * at each ice time step (every nn_fsbc time step):
310	!! - compute the modulus of ice-ocean relative velocity
311	!! at T-point (C-grid) or I-point (B-grid)
312	!! tmod_io = rhoco * \| U_ice-U_oce \|
313	!! - update the modulus of stress at ocean surface
314	!! taum = frld * taum + (1-frld) * tmod_io * \| U_ice-U_oce \|
315	!! * at each ocean time step (each kt):
316	!! compute linearized ice-ocean stresses as
317	!! Utau = tmod_io * \| U_ice - pU_oce \|
318	!! using instantaneous current ocean velocity (usually before)
319	!!
320	!! NB: - the averaging operator used depends on the ice dynamics grid (cp_ice_msh='I' or 'C')
321	!! - ice-ocean rotation angle only allowed in cp_ice_msh='I' case
322	!! - here we make an approximation: taum is only computed every ice time step
323	!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
324	!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximation...
325	!!
326	!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
327	!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
328	!!---------------------------------------------------------------------
329	INTEGER , INTENT(in) :: kt ! ocean time-step index
330	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
331	!
332	INTEGER :: ji, jj ! dummy loop indices
333	REAL(wp) :: zfrldu, zat_u, zu_i, zutau_ice, zu_t, zmodt ! local scalar
334	REAL(wp) :: zfrldv, zat_v, zv_i, zvtau_ice, zv_t, zmodi ! - -
335	REAL(wp) :: zsang, zumt ! - -
336	REAL(wp), POINTER, DIMENSION(:,:) :: ztio_u, ztio_v ! ocean stress below sea-ice
337	!!---------------------------------------------------------------------
338	!
339	CALL wrk_alloc( jpi, jpj, ztio_u, ztio_v )
340	!
341	SELECT CASE( cp_ice_msh )
342	! !-----------------------!
343	CASE( 'I' ) ! B-grid ice dynamics ! I-point (i.e. F-point with sea-ice indexation)
344	! !--=--------------------!
345	!
346	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
347	!
348	!$OMP PARALLEL
349	!$OMP DO schedule(static) private(jj,ji,zu_i,zv_i)
350	DO jj = 1, jpj !* modulus of ice-ocean relative velocity at I-point
351	DO ji = 1, jpi
352	zu_i = u_ice(ji,jj) - u_oce(ji,jj) ! ice-ocean relative velocity at I-point
353	zv_i = v_ice(ji,jj) - v_oce(ji,jj)
354	tmod_io(ji,jj) = SQRT( zu_i * zu_i + zv_i * zv_i ) ! modulus of this velocity (at I-point)
355	END DO
356	END DO
357	!$OMP DO schedule(static) private(jj,ji,zumt)
358	DO jj = 1, jpjm1 !* update the modulus of stress at ocean surface (T-point)
359	DO ji = 1, jpim1 ! NO vector opt.
360	! ! modulus of U_ice-U_oce at T-point
361	zumt = 0.25_wp * ( tmod_io(ji+1,jj) + tmod_io(ji ,jj ) &
362	& + tmod_io(ji,jj+1) + tmod_io(ji+1,jj+1) )
363	! ! update the modulus of stress at ocean surface
364	taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zumt * zumt
365	END DO
366	END DO
367	!$OMP END DO NOWAIT
368	!$OMP END PARALLEL
369	CALL lbc_lnk( taum, 'T', 1. )
370	!
371	!$OMP PARALLEL DO schedule(static) private(jj, ji)
372	DO jj = 1, jpj
373	DO ji = 1, jpi
374	utau_oce(ji,jj) = utau(ji,jj) !* save the air-ocean stresses at ice time-step
375	vtau_oce(ji,jj) = vtau(ji,jj)
376	END DO
377	END DO
378	!
379	ENDIF
380	!
381	! !== at each ocean time-step ==!
382	!
383	! !* ice/ocean stress WITH a ice-ocean rotation angle at I-point
384	!$OMP PARALLEL
385	!$OMP DO schedule(static) private(jj,ji,zsang,zu_i,zv_i,zmodi)
386	DO jj = 2, jpj
387	zsang = SIGN( 1._wp, gphif(1,jj) ) * sangvg ! change the cosine angle sign in the SH
388	DO ji = 2, jpi ! NO vect. opt. possible
389	! ... ice-ocean relative velocity at I-point using instantaneous surface ocean current at u- & v-pts
390	zu_i = u_ice(ji,jj) - 0.5_wp * ( pu_oce(ji-1,jj ) + pu_oce(ji-1,jj-1) ) * tmu(ji,jj)
391	zv_i = v_ice(ji,jj) - 0.5_wp * ( pv_oce(ji ,jj-1) + pv_oce(ji-1,jj-1) ) * tmu(ji,jj)
392	! ... components of stress with a ice-ocean rotation angle
393	zmodi = rhoco * tmod_io(ji,jj)
394	ztio_u(ji,jj) = zmodi * ( cangvg * zu_i - zsang * zv_i )
395	ztio_v(ji,jj) = zmodi * ( cangvg * zv_i + zsang * zu_i )
396	END DO
397	END DO
398	! !* surface ocean stresses at u- and v-points
399	!$OMP DO schedule(static) private(jj,ji,zutau_ice,zvtau_ice,zfrldu,zfrldv)
400	DO jj = 2, jpjm1
401	DO ji = 2, jpim1 ! NO vector opt.
402	! ! ice-ocean stress at U and V-points (from I-point values)
403	zutau_ice = 0.5_wp * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
404	zvtau_ice = 0.5_wp * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
405	! ! open-ocean (lead) fraction at U- & V-points (from T-point values)
406	zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
407	zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
408	! ! update the surface ocean stress (ice-cover wheighted)
409	utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
410	vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
411	END DO
412	END DO
413	!$OMP END DO NOWAIT
414	!$OMP END PARALLEL
415	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
416	!
417	!
418	! !-----------------------!
419	CASE( 'C' ) ! C-grid ice dynamics ! U & V-points (same as in the ocean)
420	! !--=--------------------!
421	!
422	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
423	!
424	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_t,zv_t,zmodt)
425	DO jj = 2, jpjm1 !* modulus of the ice-ocean velocity at T-point
426	DO ji = fs_2, fs_jpim1
427	zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! 2*(U_ice-U_oce) at T-point
428	zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
429	zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t ) ! \|U_ice-U_oce\|^2
430	! ! update the modulus of stress at ocean surface
431	taum (ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zmodt
432	tmod_io(ji,jj) = SQRT( zmodt ) * rhoco ! rhoco*\|Uice-Uoce\|
433	END DO
434	END DO
435	CALL lbc_lnk( taum, 'T', 1. ) ; CALL lbc_lnk( tmod_io, 'T', 1. )
436	!
437	!$OMP PARALLEL DO schedule(static) private(jj, ji)
438	DO jj = 1, jpj
439	DO ji = 1, jpi
440	utau_oce(ji,jj) = utau(ji,jj) !* save the air-ocean stresses at ice time-step
441	vtau_oce(ji,jj) = vtau(ji,jj)
442	END DO
443	END DO
444	!
445	ENDIF
446	!
447	! !== at each ocean time-step ==!
448	!
449	!$OMP PARALLEL DO schedule(static) private(jj,ji,zutau_ice,zvtau_ice,zfrldu,zfrldv)
450	DO jj = 2, jpjm1 !* ice stress over ocean WITHOUT a ice-ocean rotation angle
451	DO ji = fs_2, fs_jpim1
452	! ! ocean area at u- & v-points
453	zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
454	zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
455	! ! quadratic drag formulation without rotation
456	! ! using instantaneous surface ocean current
457	zutau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
458	zvtau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
459	! ! update the surface ocean stress (ice-cover wheighted)
460	utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
461	vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
462	END DO
463	END DO
464	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
465	!
466	END SELECT
467
468	IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, &
469	& tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask )
470	!
471	CALL wrk_dealloc( jpi, jpj, ztio_u, ztio_v )
472	!
473	END SUBROUTINE lim_sbc_tau_2
474
475
476	SUBROUTINE lim_sbc_init_2
477	!!-------------------------------------------------------------------
478	!! * ROUTINE lim_sbc_init *
479	!!
480	!! ** Purpose : Preparation of the file ice_evolu for the output of
481	!! the temporal evolution of key variables
482	!!
483	!! ** input : Namelist namicedia
484	!!-------------------------------------------------------------------
485	INTEGER :: jk,jj,ji ! local integer
486	!!-------------------------------------------------------------------
487	!
488	IF(lwp) WRITE(numout,*)
489	IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition'
490	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ '
491
492	! ! allocate lim_sbc arrays
493	IF( lim_sbc_alloc_2() /= 0 ) CALL ctl_stop( 'STOP', 'lim_sbc_flx_2 : unable to allocate arrays' )
494	!
495	r1_rdtice = 1._wp / rdt_ice
496	!
497	!$OMP PARALLEL DO schedule(static) private(jj, ji)
498	DO jj = 1, jpj
499	DO ji = 1, jpi
500	soce_0(ji,jj) = soce ! constant SSS and ice salinity used in levitating sea-ice case
501	sice_0(ji,jj) = sice
502	END DO
503	END DO
504	!
505	IF( cp_cfg == "orca" ) THEN ! decrease ocean & ice reference salinities in the Baltic sea
506	WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. &
507	& 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
508	soce_0(:,:) = 4._wp
509	sice_0(:,:) = 2._wp
510	END WHERE
511	ENDIF
512	! ! embedded sea ice
513	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
514	!$OMP PARALLEL DO schedule(static) private(jj, ji)
515	DO jj = 1, jpj
516	DO ji = 1, jpi
517	snwice_mass (ji,jj) = tms(ji,jj) * ( rhosn * hsnif(ji,jj) + rhoic * hicif(ji,jj) ) * ( 1.0 - frld(ji,jj) )
518	snwice_mass_b(ji,jj) = snwice_mass(ji,jj)
519	END DO
520	END DO
521	ELSE
522	!$OMP PARALLEL DO schedule(static) private(jj, ji)
523	DO jj = 1, jpj
524	DO ji = 1, jpi
525	snwice_mass (ji,jj) = 0.e0 ! no mass exchanges
526	snwice_mass_b(ji,jj) = 0.e0 ! no mass exchanges
527	snwice_fmass (ji,jj) = 0.e0 ! no mass exchanges
528	END DO
529	END DO
530	ENDIF
531	IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart :
532	& .NOT.ln_rstart ) THEN ! deplete the initial ssh below sea-ice area
533	!$OMP PARALLEL DO schedule(static) private(jj, ji)
534	DO jj = 1, jpj
535	DO ji = 1, jpi
536	sshn(ji,jj) = sshn(ji,jj) - snwice_mass(ji,jj) * r1_rau0
537	sshb(ji,jj) = sshb(ji,jj) - snwice_mass(ji,jj) * r1_rau0
538	END DO
539	END DO
540	!!gm I really don't like this staff here... Find a way to put that elsewhere or differently
541	!!gm
542	IF( .NOT.ln_linssh ) THEN
543
544	!$OMP PARALLEL
545	!$OMP DO schedule(static) private(jk, jj, ji)
546	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
547	DO jj = 1, jpj
548	DO ji = 1, jpi
549	e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk)( 1._wp + sshn(ji,jj)tmask(ji,jj,1)/(ht_0(ji,jj) + 1.0 - tmask(ji,jj,1)) )
550	e3t_b(ji,jj,jk) = e3t_0(ji,jj,jk)( 1._wp + sshb(ji,jj)tmask(ji,jj,1)/(ht_0(ji,jj) + 1.0 - tmask(ji,jj,1)) )
551	END DO
552	END DO
553	END DO
554	!$OMP DO schedule(static) private(jk, jj, ji)
555	DO jk = 1, jpk
556	DO jj = 1, jpj
557	DO ji = 1, jpi
558	e3t_a(ji,jj,jk) = e3t_b(ji,jj,jk)
559	END DO
560	END DO
561	END DO
562	!$OMP END PARALLEL
563	! Reconstruction of all vertical scale factors at now and before time steps
564	! ! Horizontal scale factor interpolations
565	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
566	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
567	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
568	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
569	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
570	! ! Vertical scale factor interpolations
571	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
572	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
573	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
574	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
575	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
576	! ! t- and w- points depth
577	!$OMP PARALLEL
578	!$OMP DO schedule(static) private(jj, ji)
579	DO jj = 1, jpj
580	DO ji = 1, jpi
581	gdept_n(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1)
582	gdepw_n(ji,jj,1) = 0.0_wp
583	gde3w_n(ji,jj,1) = gdept_n(ji,jj,1) - sshn(ji,jj)
584	END DO
585	END DO
586	DO jk = 2, jpk
587	!$OMP DO schedule(static) private(jj,ji)
588	DO jj = 1, jpj
589	DO ji = 1, jpi
590	gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
591	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
592	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj)
593	END DO
594	END DO
595	!$OMP END DO NOWAIT
596	END DO
597	!$OMP END PARALLEL
598	ENDIF
599	!!gm end
600	ENDIF
601	!
602	END SUBROUTINE lim_sbc_init_2
603
604	#else
605	!!----------------------------------------------------------------------
606	!! Default option Empty module NO LIM 2.0 sea-ice model
607	!!----------------------------------------------------------------------
608	#endif
609
610	!!======================================================================
611	END MODULE limsbc_2

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

Download in other formats: