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 @ 7508

Last change on this file since 7508 was 7508, checked in by mocavero, 7 years ago
changes on code duplication and workshare construct
Property svn:keywords set to `Id`
File size: 32.5 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	!$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)
139	DO jj = 1, jpj
140	DO ji = 1, jpi
141	zinda = 1.0 - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) ) ) )
142	ifvt = zinda * MAX( rzero , SIGN( rone, - phicif(ji,jj) ) )
143	i1mfr = 1.0 - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) ) ) )
144	idfr = 1.0 - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
145	iflt = zinda * (1 - i1mfr) * (1 - ifvt )
146	ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
147	iadv = ( 1 - i1mfr ) * zinda
148	ifral = ( 1 - i1mfr * ( 1 - ial ) )
149	ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
150
151	!!$ attempt to explain the tricky flags set above....
152	!!$ zinda = 1.0 - AINT( pfrld(ji,jj) ) ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo)
153	!!$ i1mfr = 1.0 - AINT( frld(ji,jj) ) ! = 0. if ice-free ocean else 1. (after ice thermo)
154	!!$
155	!!$ IF( phicif(ji,jj) <= 0. ) THEN ; ifvt = zinda ! = zinda if previous thermodynamic step overmelted the ice???
156	!!$ ELSE ; ifvt = 0. !
157	!!$ ENDIF
158	!!$
159	!!$ IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN ; idfr = 0. ! = 0. if lead fraction increases due to ice thermodynamics
160	!!$ ELSE ; idfr = 1.
161	!!$ ENDIF
162	!!$
163	!!$ iflt = zinda * (1 - i1mfr) * (1 - ifvt ) ! = 1. if ice (not only snow) at previous time and ice-free ocean currently
164	!!$
165	!!$ ial = ifvt * i1mfr + ( 1 - ifvt ) * idfr
166	!!$ = i1mfr if ifvt = 1 i.e.
167	!!$ = idfr if ifvt = 0
168	!!$! snow no ice ice ice or nothing lead fraction increases
169	!!$! at previous now at previous
170	!!$! -> ice area increases ??? -> ice area decreases ???
171	!!$
172	!!$ iadv = ( 1 - i1mfr ) * zinda
173	!!$! pure ocean ice at
174	!!$! at current previous
175	!!$! -> = 1. if ice disapear between previous and current
176	!!$
177	!!$ ifral = ( 1 - i1mfr * ( 1 - ial ) )
178	!!$! ice at ???
179	!!$! current
180	!!$! -> ???
181	!!$
182	!!$ ifrdv = ( 1 - ifral * ( 1 - ial ) ) * iadv
183	!!$! ice disapear
184	!!$
185	!!$
186
187	! computation the solar flux at ocean surface
188	IF( ln_cpl ) THEN
189	zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
190	ELSE
191	zqsr = pfrld(ji,jj) * qsr(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
192	ENDIF
193	! computation the non solar heat flux at ocean surface
194	zqns = - ( 1. - thcm(ji,jj) ) * zqsr & ! part of the solar energy used in leads
195	& + iflt * ( fscmbq(ji,jj) + ffltbif(ji,jj) ) &
196	& + ifral * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice &
197	& + ifrdv * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * r1_rdtice
198
199	fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj) ! store residual heat flux (to put into the ocean at the next time-step)
200	zqhc = ( rdq_snw(ji,jj) &
201	& + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice ! heat flux due to snow ( & ice heat content,
202	! ! if ice/ocean mass exchange active)
203	qsr (ji,jj) = zqsr ! solar heat flux
204	qns (ji,jj) = zqns - fdtcn(ji,jj) + zqhc ! non solar heat flux
205	!
206	! !------------------------------------------!
207	! ! mass and salt flux at the ocean surface !
208	! !------------------------------------------!
209	!
210	! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area)
211	! ! coupled mode:
212	IF( ln_cpl ) THEN
213	zemp = + emp_tot(ji,jj) & ! net mass flux over the grid cell (ice+ocean area)
214	& - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) ) ! minus the mass flux intercepted by sea-ice
215	ELSE
216	! ! forced mode:
217	zemp = + emp(ji,jj) * frld(ji,jj) & ! mass flux over open ocean fraction
218	& - tprecip(ji,jj) * ( 1. - frld(ji,jj) ) & ! liquid precip. over ice reaches directly the ocean
219	& + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) ) ! snow is intercepted by sea-ice (previous frld)
220	ENDIF
221	!
222	! mass flux at the ocean/ice interface (sea ice fraction)
223	zemp_snw = rdm_snw(ji,jj) * r1_rdtice ! snow melting = pure water that enters the ocean
224	zfmm = rdm_ice(ji,jj) * r1_rdtice ! Freezing minus Melting (F-M)
225
226	fmmflx(ji,jj) = zfmm ! F/M mass flux save at least for biogeochemical model
227
228	! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s]
229	zfsalt = - sice_0(ji,jj) * zfmm ! F-M salt exchange
230	zcd = soce_0(ji,jj) * zfmm ! concentration/dilution term due to F-M
231	!
232	! salt flux only : add concentration dilution term in salt flux and no F-M term in volume flux
233	! salt and mass fluxes : non concentration dilution term in salt flux and add F-M term in volume flux
234	sfx (ji,jj) = zfsalt + zswitch * zcd ! salt flux (+ C/D if no ice/ocean mass exchange)
235	emp (ji,jj) = zemp + zemp_snw + ( 1.- zswitch) * zfmm ! mass flux (+ F/M mass flux if ice/ocean mass exchange)
236	!
237	END DO
238	END DO
239	!$OMP END DO NOWAIT
240	!$OMP END PARALLEL
241	! !------------------------------------------!
242	! ! mass of snow and ice per unit area !
243	! !------------------------------------------!
244	IF( nn_ice_embd /= 0 ) THEN ! embedded sea-ice (mass required)
245	!$OMP PARALLEL DO schedule(static) private(jj, ji)
246	DO jj = 1, jpj
247	DO ji = 1, jpi
248	snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step
249	! ! new mass per unit area
250	snwice_mass (ji,jj) = tms(ji,jj) * ( rhosn * hsnif(ji,jj) + rhoic * hicif(ji,jj) ) * ( 1.0 - frld(ji,jj) )
251	! ! time evolution of snow+ice mass
252	snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) / rdt_ice
253	END DO
254	END DO
255	ENDIF
256
257	IF( iom_use('hflx_ice_cea' ) ) CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )
258	IF( iom_use('qns_io_cea' ) ) CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )
259	IF( iom_use('qsr_io_cea' ) ) CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
260
261	IF( iom_use('isnwmlt_cea' ) ) CALL iom_put( 'isnwmlt_cea' , rdm_snw(:,:) * r1_rdtice )
262	IF( iom_use('fsal_virt_cea') ) CALL iom_put( 'fsal_virt_cea', soce_0(:,:) * rdm_ice(:,:) * r1_rdtice )
263	IF( iom_use('fsal_real_cea') ) CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice )
264
265	!-----------------------------------------------!
266	! Coupling variables !
267	!-----------------------------------------------!
268
269	IF( ln_cpl) THEN
270	!$OMP PARALLEL DO schedule(static) private(jj, ji)
271	DO jj = 1, jpj
272	DO ji = 1, jpi
273	tn_ice(ji,jj,1) = sist(ji,jj) ! sea-ice surface temperature
274	ht_i(ji,jj,1) = hicif(ji,jj)
275	ht_s(ji,jj,1) = hsnif(ji,jj)
276	a_i(ji,jj,1) = fr_i(ji,jj)
277	END DO
278	END DO
279	! ! Computation of snow/ice and ocean albedo
280	CALL albedo_ice( tn_ice, ht_i, ht_s, zalbp, zalb )
281	!$OMP PARALLEL DO schedule(static) private(jj, ji)
282	DO jj = 1, jpj
283	DO ji = 1, jpi
284	alb_ice(ji,jj,1) = 0.5 * ( zalbp(ji,jj,1) + zalb (ji,jj,1) ) ! Ice albedo (mean clear and overcast skys)
285	END DO
286	END DO
287	IF( iom_use('icealb_cea' ) ) CALL iom_put( 'icealb_cea', alb_ice(:,:,1) * fr_i(:,:) ) ! ice albedo
288	ENDIF
289
290	IF(ln_ctl) THEN ! control print
291	CALL prt_ctl(tab2d_1=qsr , clinfo1=' lim_sbc: qsr : ', tab2d_2=qns , clinfo2=' qns : ')
292	CALL prt_ctl(tab2d_1=emp , clinfo1=' lim_sbc: emp : ', tab2d_2=sfx , clinfo2=' sfx : ')
293	CALL prt_ctl(tab2d_1=fr_i , clinfo1=' lim_sbc: fr_i : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice : ')
294	ENDIF
295	!
296	CALL wrk_dealloc( jpi, jpj, zqnsoce )
297	CALL wrk_dealloc( jpi, jpj, 1, zalb, zalbp )
298	!
299	END SUBROUTINE lim_sbc_flx_2
300
301
302	SUBROUTINE lim_sbc_tau_2( kt , pu_oce, pv_oce )
303	!!-------------------------------------------------------------------
304	!! * ROUTINE lim_sbc_tau_2 *
305	!!
306	!! ** Purpose : Update the ocean surface stresses due to the ice
307	!!
308	!! ** Action : * at each ice time step (every nn_fsbc time step):
309	!! - compute the modulus of ice-ocean relative velocity
310	!! at T-point (C-grid) or I-point (B-grid)
311	!! tmod_io = rhoco * \| U_ice-U_oce \|
312	!! - update the modulus of stress at ocean surface
313	!! taum = frld * taum + (1-frld) * tmod_io * \| U_ice-U_oce \|
314	!! * at each ocean time step (each kt):
315	!! compute linearized ice-ocean stresses as
316	!! Utau = tmod_io * \| U_ice - pU_oce \|
317	!! using instantaneous current ocean velocity (usually before)
318	!!
319	!! NB: - the averaging operator used depends on the ice dynamics grid (cp_ice_msh='I' or 'C')
320	!! - ice-ocean rotation angle only allowed in cp_ice_msh='I' case
321	!! - here we make an approximation: taum is only computed every ice time step
322	!! This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
323	!! to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximation...
324	!!
325	!! ** Outputs : - utau, vtau : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
326	!! - taum : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
327	!!---------------------------------------------------------------------
328	INTEGER , INTENT(in) :: kt ! ocean time-step index
329	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pu_oce, pv_oce ! surface ocean currents
330	!
331	INTEGER :: ji, jj ! dummy loop indices
332	REAL(wp) :: zfrldu, zat_u, zu_i, zutau_ice, zu_t, zmodt ! local scalar
333	REAL(wp) :: zfrldv, zat_v, zv_i, zvtau_ice, zv_t, zmodi ! - -
334	REAL(wp) :: zsang, zumt ! - -
335	REAL(wp), POINTER, DIMENSION(:,:) :: ztio_u, ztio_v ! ocean stress below sea-ice
336	!!---------------------------------------------------------------------
337	!
338	CALL wrk_alloc( jpi, jpj, ztio_u, ztio_v )
339	!
340	SELECT CASE( cp_ice_msh )
341	! !-----------------------!
342	CASE( 'I' ) ! B-grid ice dynamics ! I-point (i.e. F-point with sea-ice indexation)
343	! !--=--------------------!
344	!
345	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
346	!
347	!$OMP PARALLEL
348	!$OMP DO schedule(static) private(jj,ji,zu_i,zv_i)
349	DO jj = 1, jpj !* modulus of ice-ocean relative velocity at I-point
350	DO ji = 1, jpi
351	zu_i = u_ice(ji,jj) - u_oce(ji,jj) ! ice-ocean relative velocity at I-point
352	zv_i = v_ice(ji,jj) - v_oce(ji,jj)
353	tmod_io(ji,jj) = SQRT( zu_i * zu_i + zv_i * zv_i ) ! modulus of this velocity (at I-point)
354	END DO
355	END DO
356	!$OMP DO schedule(static) private(jj,ji,zumt)
357	DO jj = 1, jpjm1 !* update the modulus of stress at ocean surface (T-point)
358	DO ji = 1, jpim1 ! NO vector opt.
359	! ! modulus of U_ice-U_oce at T-point
360	zumt = 0.25_wp * ( tmod_io(ji+1,jj) + tmod_io(ji ,jj ) &
361	& + tmod_io(ji,jj+1) + tmod_io(ji+1,jj+1) )
362	! ! update the modulus of stress at ocean surface
363	taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zumt * zumt
364	END DO
365	END DO
366	!$OMP END DO NOWAIT
367	!$OMP END PARALLEL
368	CALL lbc_lnk( taum, 'T', 1. )
369	!
370	!$OMP PARALLEL DO schedule(static) private(jj, ji)
371	DO jj = 1, jpj
372	DO ji = 1, jpi
373	utau_oce(ji,jj) = utau(ji,jj) !* save the air-ocean stresses at ice time-step
374	vtau_oce(ji,jj) = vtau(ji,jj)
375	END DO
376	END DO
377	!
378	ENDIF
379	!
380	! !== at each ocean time-step ==!
381	!
382	! !* ice/ocean stress WITH a ice-ocean rotation angle at I-point
383	!$OMP PARALLEL
384	!$OMP DO schedule(static) private(jj,ji,zsang,zu_i,zv_i,zmodi)
385	DO jj = 2, jpj
386	zsang = SIGN( 1._wp, gphif(1,jj) ) * sangvg ! change the cosine angle sign in the SH
387	DO ji = 2, jpi ! NO vect. opt. possible
388	! ... ice-ocean relative velocity at I-point using instantaneous surface ocean current at u- & v-pts
389	zu_i = u_ice(ji,jj) - 0.5_wp * ( pu_oce(ji-1,jj ) + pu_oce(ji-1,jj-1) ) * tmu(ji,jj)
390	zv_i = v_ice(ji,jj) - 0.5_wp * ( pv_oce(ji ,jj-1) + pv_oce(ji-1,jj-1) ) * tmu(ji,jj)
391	! ... components of stress with a ice-ocean rotation angle
392	zmodi = rhoco * tmod_io(ji,jj)
393	ztio_u(ji,jj) = zmodi * ( cangvg * zu_i - zsang * zv_i )
394	ztio_v(ji,jj) = zmodi * ( cangvg * zv_i + zsang * zu_i )
395	END DO
396	END DO
397	! !* surface ocean stresses at u- and v-points
398	!$OMP DO schedule(static) private(jj,ji,zutau_ice,zvtau_ice,zfrldu,zfrldv)
399	DO jj = 2, jpjm1
400	DO ji = 2, jpim1 ! NO vector opt.
401	! ! ice-ocean stress at U and V-points (from I-point values)
402	zutau_ice = 0.5_wp * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
403	zvtau_ice = 0.5_wp * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
404	! ! open-ocean (lead) fraction at U- & V-points (from T-point values)
405	zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
406	zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
407	! ! update the surface ocean stress (ice-cover wheighted)
408	utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
409	vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
410	END DO
411	END DO
412	!$OMP END DO NOWAIT
413	!$OMP END PARALLEL
414	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
415	!
416	!
417	! !-----------------------!
418	CASE( 'C' ) ! C-grid ice dynamics ! U & V-points (same as in the ocean)
419	! !--=--------------------!
420	!
421	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN !== Ice time-step only ==! (i.e. surface module time-step)
422	!
423	!$OMP PARALLEL DO schedule(static) private(jj,ji,zu_t,zv_t,zmodt)
424	DO jj = 2, jpjm1 !* modulus of the ice-ocean velocity at T-point
425	DO ji = fs_2, fs_jpim1
426	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
427	zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
428	zmodt = 0.25_wp * ( zu_t * zu_t + zv_t * zv_t ) ! \|U_ice-U_oce\|^2
429	! ! update the modulus of stress at ocean surface
430	taum (ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zmodt
431	tmod_io(ji,jj) = SQRT( zmodt ) * rhoco ! rhoco*\|Uice-Uoce\|
432	END DO
433	END DO
434	CALL lbc_lnk( taum, 'T', 1. ) ; CALL lbc_lnk( tmod_io, 'T', 1. )
435	!
436	!$OMP PARALLEL DO schedule(static) private(jj, ji)
437	DO jj = 1, jpj
438	DO ji = 1, jpi
439	utau_oce(ji,jj) = utau(ji,jj) !* save the air-ocean stresses at ice time-step
440	vtau_oce(ji,jj) = vtau(ji,jj)
441	END DO
442	END DO
443	!
444	ENDIF
445	!
446	! !== at each ocean time-step ==!
447	!
448	!$OMP PARALLEL DO schedule(static) private(jj,ji,zutau_ice,zvtau_ice,zfrldu,zfrldv)
449	DO jj = 2, jpjm1 !* ice stress over ocean WITHOUT a ice-ocean rotation angle
450	DO ji = fs_2, fs_jpim1
451	! ! ocean area at u- & v-points
452	zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj) )
453	zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji,jj+1) )
454	! ! quadratic drag formulation without rotation
455	! ! using instantaneous surface ocean current
456	zutau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
457	zvtau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
458	! ! update the surface ocean stress (ice-cover wheighted)
459	utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
460	vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
461	END DO
462	END DO
463	CALL lbc_lnk( utau, 'U', -1. ) ; CALL lbc_lnk( vtau, 'V', -1. ) ! lateral boundary condition
464	!
465	END SELECT
466
467	IF(ln_ctl) CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau : ', mask1=umask, &
468	& tab2d_2=vtau, clinfo2=' vtau : ' , mask2=vmask )
469	!
470	CALL wrk_dealloc( jpi, jpj, ztio_u, ztio_v )
471	!
472	END SUBROUTINE lim_sbc_tau_2
473
474
475	SUBROUTINE lim_sbc_init_2
476	!!-------------------------------------------------------------------
477	!! * ROUTINE lim_sbc_init *
478	!!
479	!! ** Purpose : Preparation of the file ice_evolu for the output of
480	!! the temporal evolution of key variables
481	!!
482	!! ** input : Namelist namicedia
483	!!-------------------------------------------------------------------
484	INTEGER :: jk,jj,ji ! local integer
485	!!-------------------------------------------------------------------
486	!
487	IF(lwp) WRITE(numout,*)
488	IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition'
489	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~ '
490
491	! ! allocate lim_sbc arrays
492	IF( lim_sbc_alloc_2() /= 0 ) CALL ctl_stop( 'STOP', 'lim_sbc_flx_2 : unable to allocate arrays' )
493	!
494	r1_rdtice = 1._wp / rdt_ice
495	!
496	!$OMP PARALLEL DO schedule(static) private(jj, ji)
497	DO jj = 1, jpj
498	DO ji = 1, jpi
499	soce_0(ji,jj) = soce ! constant SSS and ice salinity used in levitating sea-ice case
500	sice_0(ji,jj) = sice
501	END DO
502	END DO
503	!
504	IF( cp_cfg == "orca" ) THEN ! decrease ocean & ice reference salinities in the Baltic sea
505	WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND. &
506	& 54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
507	soce_0(:,:) = 4._wp
508	sice_0(:,:) = 2._wp
509	END WHERE
510	ENDIF
511	! ! embedded sea ice
512	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
513	!$OMP PARALLEL DO schedule(static) private(jj, ji)
514	DO jj = 1, jpj
515	DO ji = 1, jpi
516	snwice_mass (ji,jj) = tms(ji,jj) * ( rhosn * hsnif(ji,jj) + rhoic * hicif(ji,jj) ) * ( 1.0 - frld(ji,jj) )
517	snwice_mass_b(ji,jj) = snwice_mass(ji,jj)
518	END DO
519	END DO
520	ELSE
521	!$OMP PARALLEL DO schedule(static) private(jj, ji)
522	DO jj = 1, jpj
523	DO ji = 1, jpi
524	snwice_mass (ji,jj) = 0.e0 ! no mass exchanges
525	snwice_mass_b(ji,jj) = 0.e0 ! no mass exchanges
526	snwice_fmass (ji,jj) = 0.e0 ! no mass exchanges
527	END DO
528	END DO
529	ENDIF
530	IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart :
531	& .NOT.ln_rstart ) THEN ! deplete the initial ssh below sea-ice area
532	!$OMP PARALLEL DO schedule(static) private(jj, ji)
533	DO jj = 1, jpj
534	DO ji = 1, jpi
535	sshn(ji,jj) = sshn(ji,jj) - snwice_mass(ji,jj) * r1_rau0
536	sshb(ji,jj) = sshb(ji,jj) - snwice_mass(ji,jj) * r1_rau0
537	END DO
538	END DO
539	!!gm I really don't like this staff here... Find a way to put that elsewhere or differently
540	!!gm
541	IF( .NOT.ln_linssh ) THEN
542
543	!$OMP PARALLEL
544	!$OMP DO schedule(static) private(jk)
545	do jk = 1,jpkm1 ! adjust initial vertical scale factors
546	e3t_n(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshn(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
547	e3t_b(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshb(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
548	end do
549	!$OMP DO schedule(static) private(jk, jj, ji)
550	DO jk = 1, jpk
551	DO jj = 1, jpj
552	DO ji = 1, jpi
553	e3t_a(ji,jj,jk) = e3t_b(ji,jj,jk)
554	END DO
555	END DO
556	END DO
557	!$OMP END PARALLEL
558	! Reconstruction of all vertical scale factors at now and before time steps
559	! ! Horizontal scale factor interpolations
560	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
561	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
562	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
563	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
564	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
565	! ! Vertical scale factor interpolations
566	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
567	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
568	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
569	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
570	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
571	! ! t- and w- points depth
572	!$OMP PARALLEL
573	!$OMP DO schedule(static) private(jj, ji)
574	DO jj = 1, jpj
575	DO ji = 1, jpi
576	gdept_n(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1)
577	gdepw_n(ji,jj,1) = 0.0_wp
578	gde3w_n(ji,jj,1) = gdept_n(ji,jj,1) - sshn(ji,jj)
579	END DO
580	END DO
581	DO jk = 2, jpk
582	!$OMP DO schedule(static) private(jj,ji)
583	DO jj = 1, jpj
584	DO ji = 1, jpi
585	gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
586	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
587	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj)
588	END DO
589	END DO
590	!$OMP END DO NOWAIT
591	END DO
592	!$OMP END PARALLEL
593	ENDIF
594	!!gm end
595	ENDIF
596	!
597	END SUBROUTINE lim_sbc_init_2
598
599	#else
600	!!----------------------------------------------------------------------
601	!! Default option Empty module NO LIM 2.0 sea-ice model
602	!!----------------------------------------------------------------------
603	#endif
604
605	!!======================================================================
606	END MODULE limsbc_2

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

Download in other formats: