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.

obcrad.F90 in branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

Context Navigation

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcrad.F90 @ 2797

Last change on this file since 2797 was 2797, checked in by davestorkey, 13 years ago

Delete BDY module and first implementation of new OBC module.

Initial restructuring.
Use fldread to read open boundary data.

Property svn:keywords set to Id

File size: 51.3 KB

Line
1	MODULE obcrad
2	!!=================================================================================
3	!! * MODULE obcrad *
4	!! Ocean dynamic : Phase velocities for each open boundary
5	!!=================================================================================
6	#if defined key_obc
7	!!$ !!---------------------------------------------------------------------------------
8	!!$ !! obc_rad : call the subroutine for each open boundary
9	!!$ !! obc_rad_east : compute the east phase velocities
10	!!$ !! obc_rad_west : compute the west phase velocities
11	!!$ !! obc_rad_north : compute the north phase velocities
12	!!$ !! obc_rad_south : compute the south phase velocities
13	!!$ !!---------------------------------------------------------------------------------
14	!!$ USE oce ! ocean dynamics and tracers variables
15	!!$ USE dom_oce ! ocean space and time domain variables
16	!!$ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
17	!!$ USE phycst ! physical constants
18	!!$ USE obc_oce ! ocean open boundary conditions
19	!!$ USE lib_mpp ! for mppobc
20	!!$ USE in_out_manager ! I/O units
21	!!$
22	!!$ IMPLICIT NONE
23	!!$ PRIVATE
24	!!$
25	!!$ PUBLIC obc_rad ! routine called by step.F90
26	!!$
27	!!$ INTEGER :: ji, jj, jk ! dummy loop indices
28	!!$
29	!!$ INTEGER :: & ! ... boundary space indices
30	!!$ nib = 1, & ! nib = boundary point
31	!!$ nibm = 2, & ! nibm = 1st interior point
32	!!$ nibm2 = 3, & ! nibm2 = 2nd interior point
33	!!$ ! ... boundary time indices
34	!!$ nit = 1, & ! nit = now
35	!!$ nitm = 2, & ! nitm = before
36	!!$ nitm2 = 3 ! nitm2 = before-before
37	!!$
38	!!$ !! * Substitutions
39	!!$# include "obc_vectopt_loop_substitute.h90"
40	!!$ !!---------------------------------------------------------------------------------
41	!!$ !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42	!!$ !! $Id$
43	!!$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
44	!!$ !!---------------------------------------------------------------------------------
45	!!$
46	!!$CONTAINS
47	!!$
48	!!$ SUBROUTINE obc_rad ( kt )
49	!!$ !!------------------------------------------------------------------------------
50	!!$ !! SUBROUTINE obc_rad
51	!!$ !! ********************
52	!!$ !! ** Purpose :
53	!!$ !! Perform swap of arrays to calculate radiative phase speeds at the open
54	!!$ !! boundaries and calculate those phase speeds if the open boundaries are
55	!!$ !! not fixed. In case of fixed open boundaries does nothing.
56	!!$ !!
57	!!$ !! The logical variable lp_obc_east, and/or lp_obc_west, and/or lp_obc_north,
58	!!$ !! and/or lp_obc_south allow the user to determine which boundary is an
59	!!$ !! open one (must be done in the param_obc.h90 file).
60	!!$ !!
61	!!$ !! ** Reference :
62	!!$ !! Marchesiello P., 1995, these de l'universite J. Fourier, Grenoble, France.
63	!!$ !!
64	!!$ !! History :
65	!!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90 from the
66	!!$ !! J. Molines and G. Madec version
67	!!$ !!------------------------------------------------------------------------------
68	!!$ INTEGER, INTENT( in ) :: kt
69	!!$ !!----------------------------------------------------------------------
70	!!$
71	!!$ IF( lp_obc_east .AND. .NOT.lfbceast ) CALL obc_rad_east ( kt ) ! East open boundary
72	!!$
73	!!$ IF( lp_obc_west .AND. .NOT.lfbcwest ) CALL obc_rad_west ( kt ) ! West open boundary
74	!!$
75	!!$ IF( lp_obc_north .AND. .NOT.lfbcnorth ) CALL obc_rad_north( kt ) ! North open boundary
76	!!$
77	!!$ IF( lp_obc_south .AND. .NOT.lfbcsouth ) CALL obc_rad_south( kt ) ! South open boundary
78	!!$
79	!!$ END SUBROUTINE obc_rad
80	!!$
81	!!$
82	!!$ SUBROUTINE obc_rad_east ( kt )
83	!!$ !!------------------------------------------------------------------------------
84	!!$ !! * SUBROUTINE obc_rad_east *
85	!!$ !!
86	!!$ !! ** Purpose :
87	!!$ !! Perform swap of arrays to calculate radiative phase speeds at the open
88	!!$ !! east boundary and calculate those phase speeds if this OBC is not fixed.
89	!!$ !! In case of fixed OBC, this subrountine is not called.
90	!!$ !!
91	!!$ !! History :
92	!!$ !! ! 95-03 (J.-M. Molines) Original from SPEM
93	!!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions
94	!!$ !! ! 97-12 (M. Imbard) Mpp adaptation
95	!!$ !! ! 00-06 (J.-M. Molines)
96	!!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
97	!!$ !!------------------------------------------------------------------------------
98	!!$ !! * Arguments
99	!!$ INTEGER, INTENT( in ) :: kt
100	!!$
101	!!$ !! * Local declarations
102	!!$ INTEGER :: ij
103	!!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy
104	!!$ REAL(wp) :: zucb, zucbm, zucbm2
105	!!$ !!------------------------------------------------------------------------------
106	!!$
107	!!$ ! 1. Swap arrays before calculating radiative velocities
108	!!$ ! ------------------------------------------------------
109	!!$
110	!!$ ! 1.1 zonal velocity
111	!!$ ! -------------------
112	!!$
113	!!$ IF( kt > nit000 .OR. ln_rstart ) THEN
114	!!$
115	!!$ ! ... advance in time (time filter, array swap)
116	!!$ DO jk = 1, jpkm1
117	!!$ DO jj = 1, jpj
118	!!$ uebnd(jj,jk,nib ,nitm2) = uebnd(jj,jk,nib ,nitm)*uemsk(jj,jk)
119	!!$ uebnd(jj,jk,nibm ,nitm2) = uebnd(jj,jk,nibm ,nitm)*uemsk(jj,jk)
120	!!$ uebnd(jj,jk,nibm2,nitm2) = uebnd(jj,jk,nibm2,nitm)*uemsk(jj,jk)
121	!!$ END DO
122	!!$ END DO
123	!!$ ! ... fields nitm <== nit plus time filter at the boundary
124	!!$ DO ji = fs_nie0, fs_nie1 ! Vector opt.
125	!!$ DO jk = 1, jpkm1
126	!!$ DO jj = 1, jpj
127	!!$ uebnd(jj,jk,nib ,nitm) = uebnd(jj,jk,nib, nit)*uemsk(jj,jk)
128	!!$ uebnd(jj,jk,nibm ,nitm) = uebnd(jj,jk,nibm ,nit)*uemsk(jj,jk)
129	!!$ uebnd(jj,jk,nibm2,nitm) = uebnd(jj,jk,nibm2,nit)*uemsk(jj,jk)
130	!!$ ! ... fields nit <== now (kt+1)
131	!!$ ! ... Total or baroclinic velocity at b, bm and bm2
132	!!$ zucb = un(ji,jj,jk)
133	!!$ zucbm = un(ji-1,jj,jk)
134	!!$ zucbm2 = un(ji-2,jj,jk)
135	!!$ uebnd(jj,jk,nib ,nit) = zucb *uemsk(jj,jk)
136	!!$ uebnd(jj,jk,nibm ,nit) = zucbm *uemsk(jj,jk)
137	!!$ uebnd(jj,jk,nibm2,nit) = zucbm2 *uemsk(jj,jk)
138	!!$ END DO
139	!!$ END DO
140	!!$ END DO
141	!!$ IF( lk_mpp ) CALL mppobc(uebnd,jpjed,jpjef,jpieob,jpk33,2,jpj, numout )
142	!!$
143	!!$ ! ... extremeties nie0, nie1
144	!!$ ij = jpjed +1 - njmpp
145	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
146	!!$ DO jk = 1,jpkm1
147	!!$ uebnd(ij,jk,nibm,nitm) = uebnd(ij+1 ,jk,nibm,nitm)
148	!!$ END DO
149	!!$ END IF
150	!!$ ij = jpjef +1 - njmpp
151	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
152	!!$ DO jk = 1,jpkm1
153	!!$ uebnd(ij,jk,nibm,nitm) = uebnd(ij-1 ,jk,nibm,nitm)
154	!!$ END DO
155	!!$ END IF
156	!!$
157	!!$ ! 1.2 tangential velocity
158	!!$ ! -----------------------
159	!!$
160	!!$ ! ... advance in time (time filter, array swap)
161	!!$ DO jk = 1, jpkm1
162	!!$ DO jj = 1, jpj
163	!!$ ! ... fields nitm2 <== nitm
164	!!$ vebnd(jj,jk,nib ,nitm2) = vebnd(jj,jk,nib ,nitm)*vemsk(jj,jk)
165	!!$ vebnd(jj,jk,nibm ,nitm2) = vebnd(jj,jk,nibm ,nitm)*vemsk(jj,jk)
166	!!$ vebnd(jj,jk,nibm2,nitm2) = vebnd(jj,jk,nibm2,nitm)*vemsk(jj,jk)
167	!!$ END DO
168	!!$ END DO
169	!!$
170	!!$ DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt.
171	!!$ DO jk = 1, jpkm1
172	!!$ DO jj = 1, jpj
173	!!$ vebnd(jj,jk,nib ,nitm) = vebnd(jj,jk,nib, nit)*vemsk(jj,jk)
174	!!$ vebnd(jj,jk,nibm ,nitm) = vebnd(jj,jk,nibm ,nit)*vemsk(jj,jk)
175	!!$ vebnd(jj,jk,nibm2,nitm) = vebnd(jj,jk,nibm2,nit)*vemsk(jj,jk)
176	!!$ ! ... fields nit <== now (kt+1)
177	!!$ vebnd(jj,jk,nib ,nit) = vn(ji ,jj,jk)*vemsk(jj,jk)
178	!!$ vebnd(jj,jk,nibm ,nit) = vn(ji-1,jj,jk)*vemsk(jj,jk)
179	!!$ vebnd(jj,jk,nibm2,nit) = vn(ji-2,jj,jk)*vemsk(jj,jk)
180	!!$ END DO
181	!!$ END DO
182	!!$ END DO
183	!!$ IF( lk_mpp ) CALL mppobc(vebnd,jpjed,jpjef,jpieob+1,jpk33,2,jpj, numout )
184	!!$
185	!!$ !... extremeties nie0, nie1
186	!!$ ij = jpjed +1 - njmpp
187	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
188	!!$ DO jk = 1,jpkm1
189	!!$ vebnd(ij,jk,nibm,nitm) = vebnd(ij+1 ,jk,nibm,nitm)
190	!!$ END DO
191	!!$ END IF
192	!!$ ij = jpjef +1 - njmpp
193	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
194	!!$ DO jk = 1,jpkm1
195	!!$ vebnd(ij,jk,nibm,nitm) = vebnd(ij-1 ,jk,nibm,nitm)
196	!!$ END DO
197	!!$ END IF
198	!!$
199	!!$ ! 1.3 Temperature and salinity
200	!!$ ! ----------------------------
201	!!$
202	!!$ ! ... advance in time (time filter, array swap)
203	!!$ DO jk = 1, jpkm1
204	!!$ DO jj = 1, jpj
205	!!$ ! ... fields nitm <== nit plus time filter at the boundary
206	!!$ tebnd(jj,jk,nib,nitm) = tebnd(jj,jk,nib,nit)*temsk(jj,jk)
207	!!$ sebnd(jj,jk,nib,nitm) = sebnd(jj,jk,nib,nit)*temsk(jj,jk)
208	!!$ END DO
209	!!$ END DO
210	!!$
211	!!$ DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt.
212	!!$ DO jk = 1, jpkm1
213	!!$ DO jj = 1, jpj
214	!!$ tebnd(jj,jk,nibm,nitm) = tebnd(jj,jk,nibm,nit)*temsk(jj,jk)
215	!!$ sebnd(jj,jk,nibm,nitm) = sebnd(jj,jk,nibm,nit)*temsk(jj,jk)
216	!!$ ! ... fields nit <== now (kt+1)
217	!!$ tebnd(jj,jk,nib ,nit) = tn(ji ,jj,jk)*temsk(jj,jk)
218	!!$ tebnd(jj,jk,nibm ,nit) = tn(ji-1,jj,jk)*temsk(jj,jk)
219	!!$ sebnd(jj,jk,nib ,nit) = sn(ji ,jj,jk)*temsk(jj,jk)
220	!!$ sebnd(jj,jk,nibm ,nit) = sn(ji-1,jj,jk)*temsk(jj,jk)
221	!!$ END DO
222	!!$ END DO
223	!!$ END DO
224	!!$ IF( lk_mpp ) CALL mppobc(tebnd,jpjed,jpjef,jpieob+1,jpk22,2,jpj, numout )
225	!!$ IF( lk_mpp ) CALL mppobc(sebnd,jpjed,jpjef,jpieob+1,jpk22,2,jpj, numout )
226	!!$
227	!!$ ! ... extremeties nie0, nie1
228	!!$ ij = jpjed +1 - njmpp
229	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
230	!!$ DO jk = 1,jpkm1
231	!!$ tebnd(ij,jk,nibm,nitm) = tebnd(ij+1 ,jk,nibm,nitm)
232	!!$ sebnd(ij,jk,nibm,nitm) = sebnd(ij+1 ,jk,nibm,nitm)
233	!!$ END DO
234	!!$ END IF
235	!!$ ij = jpjef +1 - njmpp
236	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
237	!!$ DO jk = 1,jpkm1
238	!!$ tebnd(ij,jk,nibm,nitm) = tebnd(ij-1 ,jk,nibm,nitm)
239	!!$ sebnd(ij,jk,nibm,nitm) = sebnd(ij-1 ,jk,nibm,nitm)
240	!!$ END DO
241	!!$ END IF
242	!!$
243	!!$ END IF ! End of array swap
244	!!$
245	!!$ ! 2 - Calculation of radiation velocities
246	!!$ ! ---------------------------------------
247	!!$
248	!!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN
249	!!$
250	!!$ ! 2.1 Calculate the normal velocity U based on phase velocity u_cxebnd
251	!!$ ! ---------------------------------------------------------------------
252	!!$ !
253	!!$ ! nibm2 nibm nib
254	!!$ ! \| nibm \| nib \|///
255	!!$ ! \| \| \| \| \|///
256	!!$ ! jj-line --f----v----f----v----f---
257	!!$ ! \| \| \| \| \|///
258	!!$ ! \| \| \|///
259	!!$ ! jj-line u T u T u///
260	!!$ ! \| \| \|///
261	!!$ ! \| \| \| \| \|///
262	!!$ ! jpieob-2 jpieob-1 jpieob
263	!!$ ! \| \|
264	!!$ ! jpieob-1 jpieob
265	!!$ !
266	!!$ ! ... (jpjedp1, jpjefm1),jpieob
267	!!$ DO ji = fs_nie0, fs_nie1 ! Vector opt.
268	!!$ DO jk = 1, jpkm1
269	!!$ DO jj = 2, jpjm1
270	!!$ ! ... 2* gradi(u) (T-point i=nibm, time mean)
271	!!$ z2dx = ( uebnd(jj,jk,nibm ,nit) + uebnd(jj,jk,nibm ,nitm2) &
272	!!$ - 2.*uebnd(jj,jk,nibm2,nitm) ) / e1t(ji-1,jj)
273	!!$ ! ... 2* gradj(u) (u-point i=nibm, time nitm)
274	!!$ z2dy = ( uebnd(jj+1,jk,nibm,nitm) - uebnd(jj-1,jk,nibm,nitm) ) / e2u(ji-1,jj)
275	!!$ ! ... square of the norm of grad(u)
276	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
277	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
278	!!$ zdt = uebnd(jj,jk,nibm,nitm2) - uebnd(jj,jk,nibm,nit)
279	!!$ ! ... i-phase speed ratio (bounded by 1)
280	!!$ IF( z4nor2 == 0. ) THEN
281	!!$ z4nor2=.00001
282	!!$ END IF
283	!!$ z05cx = zdt * z2dx / z4nor2
284	!!$ u_cxebnd(jj,jk) = z05cx*uemsk(jj,jk)
285	!!$ END DO
286	!!$ END DO
287	!!$ END DO
288	!!$
289	!!$ ! 2.2 Calculate the tangential velocity based on phase velocity v_cxebnd
290	!!$ ! -----------------------------------------------------------------------
291	!!$ !
292	!!$ ! nibm2 nibm nib
293	!!$ ! \| nibm \| nib///\|///
294	!!$ ! \| \| \| \|////\|///
295	!!$ ! jj-line --v----f----v----f----v---
296	!!$ ! \| \| \| \|////\|///
297	!!$ ! \| \| \| \|////\|///
298	!!$ ! \| jpieob-1\| jpieob /\|///
299	!!$ ! \| \| \|
300	!!$ ! jpieob-1 jpieob jpieob+1
301	!!$ !
302	!!$ ! ... (jpjedp1, jpjefm1), jpieob+1
303	!!$ DO ji = fs_nie0+1, fs_nie1+1 ! Vector opt.
304	!!$ DO jk = 1, jpkm1
305	!!$ DO jj = 2, jpjm1
306	!!$ ! ... 2* i-gradient of v (f-point i=nibm, time mean)
307	!!$ z2dx = ( vebnd(jj,jk,nibm ,nit) + vebnd(jj,jk,nibm ,nitm2) &
308	!!$ - 2.*vebnd(jj,jk,nibm2,nitm) ) / e1f(ji-2,jj)
309	!!$ ! ... 2* j-gradient of v (v-point i=nibm, time nitm)
310	!!$ z2dy = ( vebnd(jj+1,jk,nibm,nitm) - vebnd(jj-1,jk,nibm,nitm) ) / e2v(ji-1,jj)
311	!!$ ! ... square of the norm of grad(v)
312	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
313	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
314	!!$ zdt = vebnd(jj,jk,nibm,nitm2) - vebnd(jj,jk,nibm,nit)
315	!!$ ! ... i-phase speed ratio (bounded by 1) and save the unbounded phase
316	!!$ ! velocity ratio no divided by e1f for the tracer radiation
317	!!$ IF( z4nor2 == 0. ) THEN
318	!!$ z4nor2=.000001
319	!!$ END IF
320	!!$ z05cx = zdt * z2dx / z4nor2
321	!!$ v_cxebnd(jj,jk) = z05cx*vemsk(jj,jk)
322	!!$ END DO
323	!!$ END DO
324	!!$ END DO
325	!!$ IF( lk_mpp ) CALL mppobc(v_cxebnd,jpjed,jpjef,jpieob+1,jpk,2,jpj, numout )
326	!!$
327	!!$ ! ... extremeties nie0, nie1
328	!!$ ij = jpjed +1 - njmpp
329	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
330	!!$ DO jk = 1,jpkm1
331	!!$ v_cxebnd(ij,jk) = v_cxebnd(ij+1 ,jk)
332	!!$ END DO
333	!!$ END IF
334	!!$ ij = jpjef +1 - njmpp
335	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
336	!!$ DO jk = 1,jpkm1
337	!!$ v_cxebnd(ij,jk) = v_cxebnd(ij-1 ,jk)
338	!!$ END DO
339	!!$ END IF
340	!!$
341	!!$ END IF
342	!!$
343	!!$ END SUBROUTINE obc_rad_east
344	!!$
345	!!$
346	!!$ SUBROUTINE obc_rad_west ( kt )
347	!!$ !!------------------------------------------------------------------------------
348	!!$ !! * SUBROUTINE obc_rad_west *
349	!!$ !!
350	!!$ !! ** Purpose :
351	!!$ !! Perform swap of arrays to calculate radiative phase speeds at the open
352	!!$ !! west boundary and calculate those phase speeds if this OBC is not fixed.
353	!!$ !! In case of fixed OBC, this subrountine is not called.
354	!!$ !!
355	!!$ !! History :
356	!!$ !! ! 95-03 (J.-M. Molines) Original from SPEM
357	!!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions
358	!!$ !! ! 97-12 (M. Imbard) Mpp adaptation
359	!!$ !! ! 00-06 (J.-M. Molines)
360	!!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
361	!!$ !!------------------------------------------------------------------------------
362	!!$ !! * Arguments
363	!!$ INTEGER, INTENT( in ) :: kt
364	!!$
365	!!$ !! * Local declarations
366	!!$ INTEGER :: ij
367	!!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy
368	!!$ REAL(wp) :: zucb, zucbm, zucbm2
369	!!$ !!------------------------------------------------------------------------------
370	!!$
371	!!$ ! 1. Swap arrays before calculating radiative velocities
372	!!$ ! ------------------------------------------------------
373	!!$
374	!!$ ! 1.1 zonal velocity
375	!!$ ! -------------------
376	!!$
377	!!$ IF( kt > nit000 .OR. ln_rstart ) THEN
378	!!$
379	!!$ ! ... advance in time (time filter, array swap)
380	!!$ DO jk = 1, jpkm1
381	!!$ DO jj = 1, jpj
382	!!$ uwbnd(jj,jk,nib ,nitm2) = uwbnd(jj,jk,nib ,nitm)*uwmsk(jj,jk)
383	!!$ uwbnd(jj,jk,nibm ,nitm2) = uwbnd(jj,jk,nibm ,nitm)*uwmsk(jj,jk)
384	!!$ uwbnd(jj,jk,nibm2,nitm2) = uwbnd(jj,jk,nibm2,nitm)*uwmsk(jj,jk)
385	!!$ END DO
386	!!$ END DO
387	!!$
388	!!$ ! ... fields nitm <== nit plus time filter at the boundary
389	!!$ DO ji = fs_niw0, fs_niw1 ! Vector opt.
390	!!$ DO jk = 1, jpkm1
391	!!$ DO jj = 1, jpj
392	!!$ uwbnd(jj,jk,nib ,nitm) = uwbnd(jj,jk,nib ,nit)*uwmsk(jj,jk)
393	!!$ uwbnd(jj,jk,nibm ,nitm) = uwbnd(jj,jk,nibm ,nit)*uwmsk(jj,jk)
394	!!$ uwbnd(jj,jk,nibm2,nitm) = uwbnd(jj,jk,nibm2,nit)*uwmsk(jj,jk)
395	!!$ ! ... total or baroclinic velocity at b, bm and bm2
396	!!$ zucb = un (ji,jj,jk)
397	!!$ zucbm = un (ji+1,jj,jk)
398	!!$ zucbm2 = un (ji+2,jj,jk)
399	!!$
400	!!$ ! ... fields nit <== now (kt+1)
401	!!$ uwbnd(jj,jk,nib ,nit) = zucb *uwmsk(jj,jk)
402	!!$ uwbnd(jj,jk,nibm ,nit) = zucbm *uwmsk(jj,jk)
403	!!$ uwbnd(jj,jk,nibm2,nit) = zucbm2*uwmsk(jj,jk)
404	!!$ END DO
405	!!$ END DO
406	!!$ END DO
407	!!$ IF( lk_mpp ) CALL mppobc(uwbnd,jpjwd,jpjwf,jpiwob,jpk33,2,jpj, numout )
408	!!$
409	!!$ ! ... extremeties niw0, niw1
410	!!$ ij = jpjwd +1 - njmpp
411	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
412	!!$ DO jk = 1,jpkm1
413	!!$ uwbnd(ij,jk,nibm,nitm) = uwbnd(ij+1 ,jk,nibm,nitm)
414	!!$ END DO
415	!!$ END IF
416	!!$ ij = jpjwf +1 - njmpp
417	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
418	!!$ DO jk = 1,jpkm1
419	!!$ uwbnd(ij,jk,nibm,nitm) = uwbnd(ij-1 ,jk,nibm,nitm)
420	!!$ END DO
421	!!$ END IF
422	!!$
423	!!$ ! 1.2 tangential velocity
424	!!$ ! -----------------------
425	!!$
426	!!$ ! ... advance in time (time filter, array swap)
427	!!$ DO jk = 1, jpkm1
428	!!$ DO jj = 1, jpj
429	!!$ ! ... fields nitm2 <== nitm
430	!!$ vwbnd(jj,jk,nib ,nitm2) = vwbnd(jj,jk,nib ,nitm)*vwmsk(jj,jk)
431	!!$ vwbnd(jj,jk,nibm ,nitm2) = vwbnd(jj,jk,nibm ,nitm)*vwmsk(jj,jk)
432	!!$ vwbnd(jj,jk,nibm2,nitm2) = vwbnd(jj,jk,nibm2,nitm)*vwmsk(jj,jk)
433	!!$ END DO
434	!!$ END DO
435	!!$
436	!!$ DO ji = fs_niw0, fs_niw1 ! Vector opt.
437	!!$ DO jk = 1, jpkm1
438	!!$ DO jj = 1, jpj
439	!!$ vwbnd(jj,jk,nib ,nitm) = vwbnd(jj,jk,nib, nit)*vwmsk(jj,jk)
440	!!$ vwbnd(jj,jk,nibm ,nitm) = vwbnd(jj,jk,nibm ,nit)*vwmsk(jj,jk)
441	!!$ vwbnd(jj,jk,nibm2,nitm) = vwbnd(jj,jk,nibm2,nit)*vwmsk(jj,jk)
442	!!$ ! ... fields nit <== now (kt+1)
443	!!$ vwbnd(jj,jk,nib ,nit) = vn(ji ,jj,jk)*vwmsk(jj,jk)
444	!!$ vwbnd(jj,jk,nibm ,nit) = vn(ji+1,jj,jk)*vwmsk(jj,jk)
445	!!$ vwbnd(jj,jk,nibm2,nit) = vn(ji+2,jj,jk)*vwmsk(jj,jk)
446	!!$ END DO
447	!!$ END DO
448	!!$ END DO
449	!!$ IF( lk_mpp ) CALL mppobc(vwbnd,jpjwd,jpjwf,jpiwob,jpk33,2,jpj, numout )
450	!!$
451	!!$ ! ... extremeties niw0, niw1
452	!!$ ij = jpjwd +1 - njmpp
453	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
454	!!$ DO jk = 1,jpkm1
455	!!$ vwbnd(ij,jk,nibm,nitm) = vwbnd(ij+1 ,jk,nibm,nitm)
456	!!$ END DO
457	!!$ END IF
458	!!$ ij = jpjwf +1 - njmpp
459	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
460	!!$ DO jk = 1,jpkm1
461	!!$ vwbnd(ij,jk,nibm,nitm) = vwbnd(ij-1 ,jk,nibm,nitm)
462	!!$ END DO
463	!!$ END IF
464	!!$
465	!!$ ! 1.3 Temperature and salinity
466	!!$ ! ----------------------------
467	!!$
468	!!$ ! ... advance in time (time filter, array swap)
469	!!$ DO jk = 1, jpkm1
470	!!$ DO jj = 1, jpj
471	!!$ ! ... fields nitm <== nit plus time filter at the boundary
472	!!$ twbnd(jj,jk,nib,nitm) = twbnd(jj,jk,nib,nit)*twmsk(jj,jk)
473	!!$ swbnd(jj,jk,nib,nitm) = swbnd(jj,jk,nib,nit)*twmsk(jj,jk)
474	!!$ END DO
475	!!$ END DO
476	!!$
477	!!$ DO ji = fs_niw0, fs_niw1 ! Vector opt.
478	!!$ DO jk = 1, jpkm1
479	!!$ DO jj = 1, jpj
480	!!$ twbnd(jj,jk,nibm ,nitm) = twbnd(jj,jk,nibm ,nit)*twmsk(jj,jk)
481	!!$ swbnd(jj,jk,nibm ,nitm) = swbnd(jj,jk,nibm ,nit)*twmsk(jj,jk)
482	!!$ ! ... fields nit <== now (kt+1)
483	!!$ twbnd(jj,jk,nib ,nit) = tn(ji ,jj,jk)*twmsk(jj,jk)
484	!!$ twbnd(jj,jk,nibm ,nit) = tn(ji+1 ,jj,jk)*twmsk(jj,jk)
485	!!$ swbnd(jj,jk,nib ,nit) = sn(ji ,jj,jk)*twmsk(jj,jk)
486	!!$ swbnd(jj,jk,nibm ,nit) = sn(ji+1 ,jj,jk)*twmsk(jj,jk)
487	!!$ END DO
488	!!$ END DO
489	!!$ END DO
490	!!$ IF( lk_mpp ) CALL mppobc(twbnd,jpjwd,jpjwf,jpiwob,jpk22,2,jpj, numout )
491	!!$ IF( lk_mpp ) CALL mppobc(swbnd,jpjwd,jpjwf,jpiwob,jpk22,2,jpj, numout )
492	!!$
493	!!$ ! ... extremeties niw0, niw1
494	!!$ ij = jpjwd +1 - njmpp
495	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
496	!!$ DO jk = 1,jpkm1
497	!!$ twbnd(ij,jk,nibm,nitm) = twbnd(ij+1 ,jk,nibm,nitm)
498	!!$ swbnd(ij,jk,nibm,nitm) = swbnd(ij+1 ,jk,nibm,nitm)
499	!!$ END DO
500	!!$ END IF
501	!!$ ij = jpjwf +1 - njmpp
502	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
503	!!$ DO jk = 1,jpkm1
504	!!$ twbnd(ij,jk,nibm,nitm) = twbnd(ij-1 ,jk,nibm,nitm)
505	!!$ swbnd(ij,jk,nibm,nitm) = swbnd(ij-1 ,jk,nibm,nitm)
506	!!$ END DO
507	!!$ END IF
508	!!$
509	!!$ END IF ! End of array swap
510	!!$
511	!!$ ! 2 - Calculation of radiation velocities
512	!!$ ! ---------------------------------------
513	!!$
514	!!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN
515	!!$
516	!!$ ! 2.1 Calculate the normal velocity U based on phase velocity u_cxwbnd
517	!!$ ! ----------------------------------------------------------------------
518	!!$ !
519	!!$ ! nib nibm nibm2
520	!!$ ! ///\| nib \| nibm \|
521	!!$ ! ///\| \| \| \| \|
522	!!$ ! ---f----v----f----v----f-- jj-line
523	!!$ ! ///\| \| \| \| \|
524	!!$ ! ///\| \| \|
525	!!$ ! ///u T u T u jj-line
526	!!$ ! ///\| \| \|
527	!!$ ! ///\| \| \| \| \|
528	!!$ ! jpiwob jpiwob+1 jpiwob+2
529	!!$ ! \| \|
530	!!$ ! jpiwob+1 jpiwob+2
531	!!$ !
532	!!$ ! ... If free surface formulation:
533	!!$ ! ... radiative conditions on the total part + relaxation toward climatology
534	!!$ ! ... (jpjwdp1, jpjwfm1), jpiwob
535	!!$ DO ji = fs_niw0, fs_niw1 ! Vector opt.
536	!!$ DO jk = 1, jpkm1
537	!!$ DO jj = 2, jpjm1
538	!!$ ! ... 2* gradi(u) (T-point i=nibm, time mean)
539	!!$ z2dx = ( - uwbnd(jj,jk,nibm ,nit) - uwbnd(jj,jk,nibm ,nitm2) &
540	!!$ + 2.*uwbnd(jj,jk,nibm2,nitm) ) / e1t(ji+2,jj)
541	!!$ ! ... 2* gradj(u) (u-point i=nibm, time nitm)
542	!!$ z2dy = ( uwbnd(jj+1,jk,nibm,nitm) - uwbnd(jj-1,jk,nibm,nitm) ) / e2u(ji+1,jj)
543	!!$ ! ... square of the norm of grad(u)
544	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
545	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
546	!!$ zdt = uwbnd(jj,jk,nibm,nitm2) - uwbnd(jj,jk,nibm,nit)
547	!!$ ! ... i-phase speed ratio (bounded by -1)
548	!!$ IF( z4nor2 == 0. ) THEN
549	!!$ z4nor2=0.00001
550	!!$ END IF
551	!!$ z05cx = zdt * z2dx / z4nor2
552	!!$ u_cxwbnd(jj,jk)=z05cx*uwmsk(jj,jk)
553	!!$ END DO
554	!!$ END DO
555	!!$ END DO
556	!!$
557	!!$ ! 2.2 Calculate the tangential velocity based on phase velocity v_cxwbnd
558	!!$ ! -----------------------------------------------------------------------
559	!!$ !
560	!!$ ! nib nibm nibm2
561	!!$ ! ///\|///nib \| nibm \| nibm2
562	!!$ ! ///\|////\| \| \| \| \| \|
563	!!$ ! ---v----f----v----f----v----f----v-- jj-line
564	!!$ ! ///\|////\| \| \| \| \| \|
565	!!$ ! ///\|////\| \| \| \| \| \|
566	!!$ ! jpiwob jpiwob+1 jpiwob+2
567	!!$ ! \| \| \|
568	!!$ ! jpiwob jpiwob+1 jpiwob+2
569	!!$ !
570	!!$ ! ... radiative condition plus Raymond-Kuo
571	!!$ ! ... (jpjwdp1, jpjwfm1),jpiwob
572	!!$ DO ji = fs_niw0, fs_niw1 ! Vector opt.
573	!!$ DO jk = 1, jpkm1
574	!!$ DO jj = 2, jpjm1
575	!!$ ! ... 2* i-gradient of v (f-point i=nibm, time mean)
576	!!$ z2dx = ( - vwbnd(jj,jk,nibm ,nit) - vwbnd(jj,jk,nibm ,nitm2) &
577	!!$ + 2.*vwbnd(jj,jk,nibm2,nitm) ) / e1f(ji+1,jj)
578	!!$ ! ... 2* j-gradient of v (v-point i=nibm, time nitm)
579	!!$ z2dy = ( vwbnd(jj+1,jk,nibm,nitm) - vwbnd(jj-1,jk,nibm,nitm) ) / e2v(ji+1,jj)
580	!!$ ! ... square of the norm of grad(v)
581	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
582	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
583	!!$ zdt = vwbnd(jj,jk,nibm,nitm2) - vwbnd(jj,jk,nibm,nit)
584	!!$ ! ... i-phase speed ratio (bounded by -1) and save the unbounded phase
585	!!$ ! velocity ratio no divided by e1f for the tracer radiation
586	!!$ IF( z4nor2 == 0) THEN
587	!!$ z4nor2=0.000001
588	!!$ endif
589	!!$ z05cx = zdt * z2dx / z4nor2
590	!!$ v_cxwbnd(jj,jk) = z05cx*vwmsk(jj,jk)
591	!!$ END DO
592	!!$ END DO
593	!!$ END DO
594	!!$ IF( lk_mpp ) CALL mppobc(v_cxwbnd,jpjwd,jpjwf,jpiwob,jpk,2,jpj, numout )
595	!!$
596	!!$ ! ... extremeties niw0, niw1
597	!!$ ij = jpjwd +1 - njmpp
598	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
599	!!$ DO jk = 1,jpkm1
600	!!$ v_cxwbnd(ij,jk) = v_cxwbnd(ij+1 ,jk)
601	!!$ END DO
602	!!$ END IF
603	!!$ ij = jpjwf +1 - njmpp
604	!!$ IF( ij >= 2 .AND. ij < jpjm1 ) THEN
605	!!$ DO jk = 1,jpkm1
606	!!$ v_cxwbnd(ij,jk) = v_cxwbnd(ij-1 ,jk)
607	!!$ END DO
608	!!$ END IF
609	!!$
610	!!$ END IF
611	!!$
612	!!$ END SUBROUTINE obc_rad_west
613	!!$
614	!!$
615	!!$ SUBROUTINE obc_rad_north ( kt )
616	!!$ !!------------------------------------------------------------------------------
617	!!$ !! * SUBROUTINE obc_rad_north *
618	!!$ !!
619	!!$ !! ** Purpose :
620	!!$ !! Perform swap of arrays to calculate radiative phase speeds at the open
621	!!$ !! north boundary and calculate those phase speeds if this OBC is not fixed.
622	!!$ !! In case of fixed OBC, this subrountine is not called.
623	!!$ !!
624	!!$ !! History :
625	!!$ !! ! 95-03 (J.-M. Molines) Original from SPEM
626	!!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions
627	!!$ !! ! 97-12 (M. Imbard) Mpp adaptation
628	!!$ !! ! 00-06 (J.-M. Molines)
629	!!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
630	!!$ !!------------------------------------------------------------------------------
631	!!$ !! * Arguments
632	!!$ INTEGER, INTENT( in ) :: kt
633	!!$
634	!!$ !! * Local declarations
635	!!$ INTEGER :: ii
636	!!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy
637	!!$ REAL(wp) :: zvcb, zvcbm, zvcbm2
638	!!$ !!------------------------------------------------------------------------------
639	!!$
640	!!$ ! 1. Swap arrays before calculating radiative velocities
641	!!$ ! ------------------------------------------------------
642	!!$
643	!!$ ! 1.1 zonal velocity
644	!!$ ! -------------------
645	!!$
646	!!$ IF( kt > nit000 .OR. ln_rstart ) THEN
647	!!$
648	!!$ ! ... advance in time (time filter, array swap)
649	!!$ DO jk = 1, jpkm1
650	!!$ DO ji = 1, jpi
651	!!$ ! ... fields nitm2 <== nitm
652	!!$ unbnd(ji,jk,nib ,nitm2) = unbnd(ji,jk,nib ,nitm)*unmsk(ji,jk)
653	!!$ unbnd(ji,jk,nibm ,nitm2) = unbnd(ji,jk,nibm ,nitm)*unmsk(ji,jk)
654	!!$ unbnd(ji,jk,nibm2,nitm2) = unbnd(ji,jk,nibm2,nitm)*unmsk(ji,jk)
655	!!$ END DO
656	!!$ END DO
657	!!$
658	!!$ DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt.
659	!!$ DO jk = 1, jpkm1
660	!!$ DO ji = 1, jpi
661	!!$ unbnd(ji,jk,nib ,nitm) = unbnd(ji,jk,nib, nit)*unmsk(ji,jk)
662	!!$ unbnd(ji,jk,nibm ,nitm) = unbnd(ji,jk,nibm ,nit)*unmsk(ji,jk)
663	!!$ unbnd(ji,jk,nibm2,nitm) = unbnd(ji,jk,nibm2,nit)*unmsk(ji,jk)
664	!!$ ! ... fields nit <== now (kt+1)
665	!!$ unbnd(ji,jk,nib ,nit) = un(ji,jj, jk)*unmsk(ji,jk)
666	!!$ unbnd(ji,jk,nibm ,nit) = un(ji,jj-1,jk)*unmsk(ji,jk)
667	!!$ unbnd(ji,jk,nibm2,nit) = un(ji,jj-2,jk)*unmsk(ji,jk)
668	!!$ END DO
669	!!$ END DO
670	!!$ END DO
671	!!$ IF( lk_mpp ) CALL mppobc(unbnd,jpind,jpinf,jpjnob+1,jpk33,1,jpi, numout )
672	!!$
673	!!$ ! ... extremeties njn0,njn1
674	!!$ ii = jpind + 1 - nimpp
675	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
676	!!$ DO jk = 1, jpkm1
677	!!$ unbnd(ii,jk,nibm,nitm) = unbnd(ii+1,jk,nibm,nitm)
678	!!$ END DO
679	!!$ END IF
680	!!$ ii = jpinf + 1 - nimpp
681	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
682	!!$ DO jk = 1, jpkm1
683	!!$ unbnd(ii,jk,nibm,nitm) = unbnd(ii-1,jk,nibm,nitm)
684	!!$ END DO
685	!!$ END IF
686	!!$
687	!!$ ! 1.2. normal velocity
688	!!$ ! --------------------
689	!!$
690	!!$ ! ... advance in time (time filter, array swap)
691	!!$ DO jk = 1, jpkm1
692	!!$ DO ji = 1, jpi
693	!!$ ! ... fields nitm2 <== nitm
694	!!$ vnbnd(ji,jk,nib ,nitm2) = vnbnd(ji,jk,nib ,nitm)*vnmsk(ji,jk)
695	!!$ vnbnd(ji,jk,nibm ,nitm2) = vnbnd(ji,jk,nibm ,nitm)*vnmsk(ji,jk)
696	!!$ vnbnd(ji,jk,nibm2,nitm2) = vnbnd(ji,jk,nibm2,nitm)*vnmsk(ji,jk)
697	!!$ END DO
698	!!$ END DO
699	!!$
700	!!$ DO jj = fs_njn0, fs_njn1 ! Vector opt.
701	!!$ DO jk = 1, jpkm1
702	!!$ DO ji = 1, jpi
703	!!$ vnbnd(ji,jk,nib ,nitm) = vnbnd(ji,jk,nib, nit)*vnmsk(ji,jk)
704	!!$ vnbnd(ji,jk,nibm ,nitm) = vnbnd(ji,jk,nibm ,nit)*vnmsk(ji,jk)
705	!!$ vnbnd(ji,jk,nibm2,nitm) = vnbnd(ji,jk,nibm2,nit)*vnmsk(ji,jk)
706	!!$ ! ... fields nit <== now (kt+1)
707	!!$ ! ... total or baroclinic velocity at b, bm and bm2
708	!!$ zvcb = vn (ji,jj,jk)
709	!!$ zvcbm = vn (ji,jj-1,jk)
710	!!$ zvcbm2 = vn (ji,jj-2,jk)
711	!!$ ! ... fields nit <== now (kt+1)
712	!!$ vnbnd(ji,jk,nib ,nit) = zvcb *vnmsk(ji,jk)
713	!!$ vnbnd(ji,jk,nibm ,nit) = zvcbm *vnmsk(ji,jk)
714	!!$ vnbnd(ji,jk,nibm2,nit) = zvcbm2*vnmsk(ji,jk)
715	!!$ END DO
716	!!$ END DO
717	!!$ END DO
718	!!$ IF( lk_mpp ) CALL mppobc(vnbnd,jpind,jpinf,jpjnob,jpk33,1,jpi, numout )
719	!!$
720	!!$ ! ... extremeties njn0,njn1
721	!!$ ii = jpind + 1 - nimpp
722	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
723	!!$ DO jk = 1, jpkm1
724	!!$ vnbnd(ii,jk,nibm,nitm) = vnbnd(ii+1,jk,nibm,nitm)
725	!!$ END DO
726	!!$ END IF
727	!!$ ii = jpinf + 1 - nimpp
728	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
729	!!$ DO jk = 1, jpkm1
730	!!$ vnbnd(ii,jk,nibm,nitm) = vnbnd(ii-1,jk,nibm,nitm)
731	!!$ END DO
732	!!$ END IF
733	!!$
734	!!$ ! 1.3 Temperature and salinity
735	!!$ ! ----------------------------
736	!!$
737	!!$ ! ... advance in time (time filter, array swap)
738	!!$ DO jk = 1, jpkm1
739	!!$ DO ji = 1, jpi
740	!!$ ! ... fields nitm <== nit plus time filter at the boundary
741	!!$ tnbnd(ji,jk,nib ,nitm) = tnbnd(ji,jk,nib,nit)*tnmsk(ji,jk)
742	!!$ snbnd(ji,jk,nib ,nitm) = snbnd(ji,jk,nib,nit)*tnmsk(ji,jk)
743	!!$ END DO
744	!!$ END DO
745	!!$
746	!!$ DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt.
747	!!$ DO jk = 1, jpkm1
748	!!$ DO ji = 1, jpi
749	!!$ tnbnd(ji,jk,nibm ,nitm) = tnbnd(ji,jk,nibm ,nit)*tnmsk(ji,jk)
750	!!$ snbnd(ji,jk,nibm ,nitm) = snbnd(ji,jk,nibm ,nit)*tnmsk(ji,jk)
751	!!$ ! ... fields nit <== now (kt+1)
752	!!$ tnbnd(ji,jk,nib ,nit) = tn(ji,jj, jk)*tnmsk(ji,jk)
753	!!$ tnbnd(ji,jk,nibm ,nit) = tn(ji,jj-1,jk)*tnmsk(ji,jk)
754	!!$ snbnd(ji,jk,nib ,nit) = sn(ji,jj, jk)*tnmsk(ji,jk)
755	!!$ snbnd(ji,jk,nibm ,nit) = sn(ji,jj-1,jk)*tnmsk(ji,jk)
756	!!$ END DO
757	!!$ END DO
758	!!$ END DO
759	!!$ IF( lk_mpp ) CALL mppobc(tnbnd,jpind,jpinf,jpjnob+1,jpk22,1,jpi, numout )
760	!!$ IF( lk_mpp ) CALL mppobc(snbnd,jpind,jpinf,jpjnob+1,jpk22,1,jpi, numout )
761	!!$
762	!!$ ! ... extremeties njn0,njn1
763	!!$ ii = jpind + 1 - nimpp
764	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
765	!!$ DO jk = 1, jpkm1
766	!!$ tnbnd(ii,jk,nibm,nitm) = tnbnd(ii+1,jk,nibm,nitm)
767	!!$ snbnd(ii,jk,nibm,nitm) = snbnd(ii+1,jk,nibm,nitm)
768	!!$ END DO
769	!!$ END IF
770	!!$ ii = jpinf + 1 - nimpp
771	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
772	!!$ DO jk = 1, jpkm1
773	!!$ tnbnd(ii,jk,nibm,nitm) = tnbnd(ii-1,jk,nibm,nitm)
774	!!$ snbnd(ii,jk,nibm,nitm) = snbnd(ii-1,jk,nibm,nitm)
775	!!$ END DO
776	!!$ END IF
777	!!$
778	!!$ END IF ! End of array swap
779	!!$
780	!!$ ! 2 - Calculation of radiation velocities
781	!!$ ! ---------------------------------------
782	!!$
783	!!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN
784	!!$
785	!!$ ! 2.1 Calculate the normal velocity based on phase velocity u_cynbnd
786	!!$ ! -------------------------------------------------------------------
787	!!$ !
788	!!$ ! ji-row
789	!!$ ! \|
790	!!$ ! nib -///u////// jpjnob + 1
791	!!$ ! /////\|//////
792	!!$ ! nib -----f----- jpjnob
793	!!$ ! \|
794	!!$ ! nibm-- u ---- jpjnob
795	!!$ ! \|
796	!!$ ! nibm -----f----- jpjnob-1
797	!!$ ! \|
798	!!$ ! nibm2-- u ---- jpjnob-1
799	!!$ ! \|
800	!!$ ! nibm2 -----f----- jpjnob-2
801	!!$ ! \|
802	!!$ ! ... radiative condition
803	!!$ ! ... jpjnob+1,(jpindp1, jpinfm1)
804	!!$ DO jj = fs_njn0+1, fs_njn1+1 ! Vector opt.
805	!!$ DO jk = 1, jpkm1
806	!!$ DO ji = 2, jpim1
807	!!$ ! ... 2* j-gradient of u (f-point i=nibm, time mean)
808	!!$ z2dx = ( unbnd(ji,jk,nibm ,nit) + unbnd(ji,jk,nibm ,nitm2) &
809	!!$ - 2.*unbnd(ji,jk,nibm2,nitm)) / e2f(ji,jj-2)
810	!!$ ! ... 2* i-gradient of u (u-point i=nibm, time nitm)
811	!!$ z2dy = ( unbnd(ji+1,jk,nibm,nitm) - unbnd(ji-1,jk,nibm,nitm) ) / e1u(ji,jj-1)
812	!!$ ! ... square of the norm of grad(v)
813	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
814	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
815	!!$ zdt = unbnd(ji,jk,nibm,nitm2) - unbnd(ji,jk,nibm,nit)
816	!!$ ! ... i-phase speed ratio (bounded by 1) and save the unbounded phase
817	!!$ ! velocity ratio no divided by e1f for the tracer radiation
818	!!$ IF( z4nor2 == 0.) THEN
819	!!$ z4nor2=.000001
820	!!$ END IF
821	!!$ z05cx = zdt * z2dx / z4nor2
822	!!$ u_cynbnd(ji,jk) = z05cx *unmsk(ji,jk)
823	!!$ END DO
824	!!$ END DO
825	!!$ END DO
826	!!$ IF( lk_mpp ) CALL mppobc(u_cynbnd,jpind,jpinf,jpjnob+1,jpk,1,jpi, numout )
827	!!$
828	!!$ ! ... extremeties njn0,njn1
829	!!$ ii = jpind + 1 - nimpp
830	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
831	!!$ DO jk = 1, jpkm1
832	!!$ u_cynbnd(ii,jk) = u_cynbnd(ii+1,jk)
833	!!$ END DO
834	!!$ END IF
835	!!$ ii = jpinf + 1 - nimpp
836	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
837	!!$ DO jk = 1, jpkm1
838	!!$ u_cynbnd(ii,jk) = u_cynbnd(ii-1,jk)
839	!!$ END DO
840	!!$ END IF
841	!!$
842	!!$ ! 2.2 Calculate the normal velocity based on phase velocity v_cynbnd
843	!!$ ! ------------------------------------------------------------------
844	!!$ !
845	!!$ ! ji-row ji-row
846	!!$ ! \|
847	!!$ ! /////\|/////////////////
848	!!$ ! nib -----f----v----f---- jpjnob
849	!!$ ! \| \|
850	!!$ ! nib - u -- T -- u ---- jpjnob
851	!!$ ! \| \|
852	!!$ ! nibm -----f----v----f---- jpjnob-1
853	!!$ ! \| \|
854	!!$ ! nibm -- u -- T -- u --- jpjnob-1
855	!!$ ! \| \|
856	!!$ ! nibm2 -----f----v----f---- jpjnob-2
857	!!$ ! \| \|
858	!!$ ! ... Free surface formulation:
859	!!$ ! ... radiative conditions on the total part + relaxation toward climatology
860	!!$ ! ... jpjnob,(jpindp1, jpinfm1)
861	!!$ DO jj = fs_njn0, fs_njn1 ! Vector opt.
862	!!$ DO jk = 1, jpkm1
863	!!$ DO ji = 2, jpim1
864	!!$ ! ... 2* gradj(v) (T-point i=nibm, time mean)
865	!!$ ii = ji -1 + nimpp
866	!!$ z2dx = ( vnbnd(ji,jk,nibm ,nit) + vnbnd(ji,jk,nibm ,nitm2) &
867	!!$ - 2.*vnbnd(ji,jk,nibm2,nitm)) / e2t(ji,jj-1)
868	!!$ ! ... 2* gradi(v) (v-point i=nibm, time nitm)
869	!!$ z2dy = ( vnbnd(ji+1,jk,nibm,nitm) - vnbnd(ji-1,jk,nibm,nitm) ) / e1v(ji,jj-1)
870	!!$ ! ... square of the norm of grad(u)
871	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
872	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
873	!!$ zdt = vnbnd(ji,jk,nibm,nitm2) - vnbnd(ji,jk,nibm,nit)
874	!!$ ! ... j-phase speed ratio (bounded by 1)
875	!!$ IF( z4nor2 == 0. ) THEN
876	!!$ z4nor2=.00001
877	!!$ END IF
878	!!$ z05cx = zdt * z2dx / z4nor2
879	!!$ v_cynbnd(ji,jk)=z05cx *vnmsk(ji,jk)
880	!!$ END DO
881	!!$ END DO
882	!!$ END DO
883	!!$
884	!!$ END IF
885	!!$
886	!!$ END SUBROUTINE obc_rad_north
887	!!$
888	!!$
889	!!$ SUBROUTINE obc_rad_south ( kt )
890	!!$ !!------------------------------------------------------------------------------
891	!!$ !! * SUBROUTINE obc_rad_south *
892	!!$ !!
893	!!$ !! ** Purpose :
894	!!$ !! Perform swap of arrays to calculate radiative phase speeds at the open
895	!!$ !! south boundary and calculate those phase speeds if this OBC is not fixed.
896	!!$ !! In case of fixed OBC, this subrountine is not called.
897	!!$ !!
898	!!$ !! History :
899	!!$ !! ! 95-03 (J.-M. Molines) Original from SPEM
900	!!$ !! ! 97-07 (G. Madec, J.-M. Molines) additions
901	!!$ !! ! 97-12 (M. Imbard) Mpp adaptation
902	!!$ !! ! 00-06 (J.-M. Molines)
903	!!$ !! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
904	!!$ !!------------------------------------------------------------------------------
905	!!$ !! * Arguments
906	!!$ INTEGER, INTENT( in ) :: kt
907	!!$
908	!!$ !! * Local declarations
909	!!$ INTEGER :: ii
910	!!$ REAL(wp) :: z05cx, zdt, z4nor2, z2dx, z2dy
911	!!$ REAL(wp) :: zvcb, zvcbm, zvcbm2
912	!!$ !!------------------------------------------------------------------------------
913	!!$
914	!!$ ! 1. Swap arrays before calculating radiative velocities
915	!!$ ! ------------------------------------------------------
916	!!$
917	!!$ ! 1.1 zonal velocity
918	!!$ ! --------------------
919	!!$
920	!!$ IF( kt > nit000 .OR. ln_rstart ) THEN
921	!!$
922	!!$ ! ... advance in time (time filter, array swap)
923	!!$ DO jk = 1, jpkm1
924	!!$ DO ji = 1, jpi
925	!!$ ! ... fields nitm2 <== nitm
926	!!$ usbnd(ji,jk,nib ,nitm2) = usbnd(ji,jk,nib ,nitm)*usmsk(ji,jk)
927	!!$ usbnd(ji,jk,nibm ,nitm2) = usbnd(ji,jk,nibm ,nitm)*usmsk(ji,jk)
928	!!$ usbnd(ji,jk,nibm2,nitm2) = usbnd(ji,jk,nibm2,nitm)*usmsk(ji,jk)
929	!!$ END DO
930	!!$ END DO
931	!!$
932	!!$ DO jj = fs_njs0, fs_njs1 ! Vector opt.
933	!!$ DO jk = 1, jpkm1
934	!!$ DO ji = 1, jpi
935	!!$ usbnd(ji,jk,nib ,nitm) = usbnd(ji,jk,nib, nit)*usmsk(ji,jk)
936	!!$ usbnd(ji,jk,nibm ,nitm) = usbnd(ji,jk,nibm ,nit)*usmsk(ji,jk)
937	!!$ usbnd(ji,jk,nibm2,nitm) = usbnd(ji,jk,nibm2,nit)*usmsk(ji,jk)
938	!!$ ! ... fields nit <== now (kt+1)
939	!!$ usbnd(ji,jk,nib ,nit) = un(ji,jj ,jk)*usmsk(ji,jk)
940	!!$ usbnd(ji,jk,nibm ,nit) = un(ji,jj+1,jk)*usmsk(ji,jk)
941	!!$ usbnd(ji,jk,nibm2,nit) = un(ji,jj+2,jk)*usmsk(ji,jk)
942	!!$ END DO
943	!!$ END DO
944	!!$ END DO
945	!!$ IF( lk_mpp ) CALL mppobc(usbnd,jpisd,jpisf,jpjsob,jpk33,1,jpi, numout )
946	!!$
947	!!$ ! ... extremeties njs0,njs1
948	!!$ ii = jpisd + 1 - nimpp
949	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
950	!!$ DO jk = 1, jpkm1
951	!!$ usbnd(ii,jk,nibm,nitm) = usbnd(ii+1,jk,nibm,nitm)
952	!!$ END DO
953	!!$ END IF
954	!!$ ii = jpisf + 1 - nimpp
955	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
956	!!$ DO jk = 1, jpkm1
957	!!$ usbnd(ii,jk,nibm,nitm) = usbnd(ii-1,jk,nibm,nitm)
958	!!$ END DO
959	!!$ END IF
960	!!$
961	!!$ ! 1.2 normal velocity
962	!!$ ! -------------------
963	!!$
964	!!$ !.. advance in time (time filter, array swap)
965	!!$ DO jk = 1, jpkm1
966	!!$ DO ji = 1, jpi
967	!!$ ! ... fields nitm2 <== nitm
968	!!$ vsbnd(ji,jk,nib ,nitm2) = vsbnd(ji,jk,nib ,nitm)*vsmsk(ji,jk)
969	!!$ vsbnd(ji,jk,nibm ,nitm2) = vsbnd(ji,jk,nibm ,nitm)*vsmsk(ji,jk)
970	!!$ END DO
971	!!$ END DO
972	!!$
973	!!$ DO jj = fs_njs0, fs_njs1 ! Vector opt.
974	!!$ DO jk = 1, jpkm1
975	!!$ DO ji = 1, jpi
976	!!$ vsbnd(ji,jk,nib ,nitm) = vsbnd(ji,jk,nib, nit)*vsmsk(ji,jk)
977	!!$ vsbnd(ji,jk,nibm ,nitm) = vsbnd(ji,jk,nibm ,nit)*vsmsk(ji,jk)
978	!!$ vsbnd(ji,jk,nibm2,nitm) = vsbnd(ji,jk,nibm2,nit)*vsmsk(ji,jk)
979	!!$ ! ... total or baroclinic velocity at b, bm and bm2
980	!!$ zvcb = vn (ji,jj,jk)
981	!!$ zvcbm = vn (ji,jj+1,jk)
982	!!$ zvcbm2 = vn (ji,jj+2,jk)
983	!!$ ! ... fields nit <== now (kt+1)
984	!!$ vsbnd(ji,jk,nib ,nit) = zvcb *vsmsk(ji,jk)
985	!!$ vsbnd(ji,jk,nibm ,nit) = zvcbm *vsmsk(ji,jk)
986	!!$ vsbnd(ji,jk,nibm2,nit) = zvcbm2 *vsmsk(ji,jk)
987	!!$ END DO
988	!!$ END DO
989	!!$ END DO
990	!!$ IF( lk_mpp ) CALL mppobc(vsbnd,jpisd,jpisf,jpjsob,jpk33,1,jpi, numout )
991	!!$
992	!!$ ! ... extremeties njs0,njs1
993	!!$ ii = jpisd + 1 - nimpp
994	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
995	!!$ DO jk = 1, jpkm1
996	!!$ vsbnd(ii,jk,nibm,nitm) = vsbnd(ii+1,jk,nibm,nitm)
997	!!$ END DO
998	!!$ END IF
999	!!$ ii = jpisf + 1 - nimpp
1000	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
1001	!!$ DO jk = 1, jpkm1
1002	!!$ vsbnd(ii,jk,nibm,nitm) = vsbnd(ii-1,jk,nibm,nitm)
1003	!!$ END DO
1004	!!$ END IF
1005	!!$
1006	!!$ ! 1.3 Temperature and salinity
1007	!!$ ! ----------------------------
1008	!!$
1009	!!$ ! ... advance in time (time filter, array swap)
1010	!!$ DO jk = 1, jpkm1
1011	!!$ DO ji = 1, jpi
1012	!!$ ! ... fields nitm <== nit plus time filter at the boundary
1013	!!$ tsbnd(ji,jk,nib,nitm) = tsbnd(ji,jk,nib,nit)*tsmsk(ji,jk)
1014	!!$ ssbnd(ji,jk,nib,nitm) = ssbnd(ji,jk,nib,nit)*tsmsk(ji,jk)
1015	!!$ END DO
1016	!!$ END DO
1017	!!$
1018	!!$ DO jj = fs_njs0, fs_njs1 ! Vector opt.
1019	!!$ DO jk = 1, jpkm1
1020	!!$ DO ji = 1, jpi
1021	!!$ tsbnd(ji,jk,nibm ,nitm) = tsbnd(ji,jk,nibm ,nit)*tsmsk(ji,jk)
1022	!!$ ssbnd(ji,jk,nibm ,nitm) = ssbnd(ji,jk,nibm ,nit)*tsmsk(ji,jk)
1023	!!$ ! ... fields nit <== now (kt+1)
1024	!!$ tsbnd(ji,jk,nib ,nit) = tn(ji,jj ,jk)*tsmsk(ji,jk)
1025	!!$ tsbnd(ji,jk,nibm ,nit) = tn(ji,jj+1 ,jk)*tsmsk(ji,jk)
1026	!!$ ssbnd(ji,jk,nib ,nit) = sn(ji,jj ,jk)*tsmsk(ji,jk)
1027	!!$ ssbnd(ji,jk,nibm ,nit) = sn(ji,jj+1 ,jk)*tsmsk(ji,jk)
1028	!!$ END DO
1029	!!$ END DO
1030	!!$ END DO
1031	!!$ IF( lk_mpp ) CALL mppobc(tsbnd,jpisd,jpisf,jpjsob,jpk22,1,jpi, numout )
1032	!!$ IF( lk_mpp ) CALL mppobc(ssbnd,jpisd,jpisf,jpjsob,jpk22,1,jpi, numout )
1033	!!$
1034	!!$ ! ... extremeties njs0,njs1
1035	!!$ ii = jpisd + 1 - nimpp
1036	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
1037	!!$ DO jk = 1, jpkm1
1038	!!$ tsbnd(ii,jk,nibm,nitm) = tsbnd(ii+1,jk,nibm,nitm)
1039	!!$ ssbnd(ii,jk,nibm,nitm) = ssbnd(ii+1,jk,nibm,nitm)
1040	!!$ END DO
1041	!!$ END IF
1042	!!$ ii = jpisf + 1 - nimpp
1043	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
1044	!!$ DO jk = 1, jpkm1
1045	!!$ tsbnd(ii,jk,nibm,nitm) = tsbnd(ii-1,jk,nibm,nitm)
1046	!!$ ssbnd(ii,jk,nibm,nitm) = ssbnd(ii-1,jk,nibm,nitm)
1047	!!$ END DO
1048	!!$ END IF
1049	!!$
1050	!!$ END IF ! End of array swap
1051	!!$
1052	!!$ ! 2 - Calculation of radiation velocities
1053	!!$ ! ---------------------------------------
1054	!!$
1055	!!$ IF( kt >= nit000 +3 .OR. ln_rstart ) THEN
1056	!!$
1057	!!$ ! 2.1 Calculate the normal velocity based on phase velocity u_cysbnd
1058	!!$ ! -------------------------------------------------------------------
1059	!!$ !
1060	!!$ ! ji-row
1061	!!$ ! \|
1062	!!$ ! nibm2 -----f----- jpjsob +2
1063	!!$ ! \|
1064	!!$ ! nibm2 -- u ----- jpjsob +2
1065	!!$ ! \|
1066	!!$ ! nibm -----f----- jpjsob +1
1067	!!$ ! \|
1068	!!$ ! nibm -- u ----- jpjsob +1
1069	!!$ ! \|
1070	!!$ ! nib -----f----- jpjsob
1071	!!$ ! /////\|//////
1072	!!$ ! nib ////u///// jpjsob
1073	!!$ !
1074	!!$ ! ... radiative condition plus Raymond-Kuo
1075	!!$ ! ... jpjsob,(jpisdp1, jpisfm1)
1076	!!$ DO jj = fs_njs0, fs_njs1 ! Vector opt.
1077	!!$ DO jk = 1, jpkm1
1078	!!$ DO ji = 2, jpim1
1079	!!$ ! ... 2* j-gradient of u (f-point i=nibm, time mean)
1080	!!$ z2dx = (- usbnd(ji,jk,nibm ,nit) - usbnd(ji,jk,nibm ,nitm2) &
1081	!!$ + 2.*usbnd(ji,jk,nibm2,nitm) ) / e2f(ji,jj+1)
1082	!!$ ! ... 2* i-gradient of u (u-point i=nibm, time nitm)
1083	!!$ z2dy = ( usbnd(ji+1,jk,nibm,nitm) - usbnd(ji-1,jk,nibm,nitm) ) / e1u(ji, jj+1)
1084	!!$ ! ... square of the norm of grad(v)
1085	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
1086	!!$ IF( z4nor2 == 0.) THEN
1087	!!$ z4nor2 = 0.000001
1088	!!$ END IF
1089	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
1090	!!$ zdt = usbnd(ji,jk,nibm,nitm2) - usbnd(ji,jk,nibm,nit)
1091	!!$ ! ... i-phase speed ratio (bounded by -1) and save the unbounded phase
1092	!!$ ! velocity ratio no divided by e1f for the tracer radiation
1093	!!$ z05cx = zdt * z2dx / z4nor2
1094	!!$ u_cysbnd(ji,jk) = z05cx*usmsk(ji,jk)
1095	!!$ END DO
1096	!!$ END DO
1097	!!$ END DO
1098	!!$ IF( lk_mpp ) CALL mppobc(u_cysbnd,jpisd,jpisf,jpjsob,jpk,1,jpi, numout )
1099	!!$
1100	!!$ ! ... extremeties njs0,njs1
1101	!!$ ii = jpisd + 1 - nimpp
1102	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
1103	!!$ DO jk = 1, jpkm1
1104	!!$ u_cysbnd(ii,jk) = u_cysbnd(ii+1,jk)
1105	!!$ END DO
1106	!!$ END IF
1107	!!$ ii = jpisf + 1 - nimpp
1108	!!$ IF( ii >= 2 .AND. ii < jpim1 ) THEN
1109	!!$ DO jk = 1, jpkm1
1110	!!$ u_cysbnd(ii,jk) = u_cysbnd(ii-1,jk)
1111	!!$ END DO
1112	!!$ END IF
1113	!!$
1114	!!$ ! 2.2 Calculate the normal velocity based on phase velocity v_cysbnd
1115	!!$ ! -------------------------------------------------------------------
1116	!!$ !
1117	!!$ ! ji-row ji-row
1118	!!$ ! \| \|
1119	!!$ ! nibm2 -----f----v----f---- jpjsob+2
1120	!!$ ! \| \|
1121	!!$ ! nibm - u -- T -- u ---- jpjsob+2
1122	!!$ ! \| \|
1123	!!$ ! nibm -----f----v----f---- jpjsob+1
1124	!!$ ! \| \|
1125	!!$ ! nib -- u -- T -- u --- jpjsob+1
1126	!!$ ! \| \|
1127	!!$ ! nib -----f----v----f---- jpjsob
1128	!!$ ! /////////////////////
1129	!!$ !
1130	!!$ ! ... Free surface formulation:
1131	!!$ ! ... radiative conditions on the total part + relaxation toward climatology
1132	!!$ ! ... jpjsob,(jpisdp1,jpisfm1)
1133	!!$ DO jj = fs_njs0, fs_njs1 ! Vector opt.
1134	!!$ DO jk = 1, jpkm1
1135	!!$ DO ji = 2, jpim1
1136	!!$ ! ... 2* gradj(v) (T-point i=nibm, time mean)
1137	!!$ z2dx = ( - vsbnd(ji,jk,nibm ,nit) - vsbnd(ji,jk,nibm ,nitm2) &
1138	!!$ + 2.*vsbnd(ji,jk,nibm2,nitm) ) / e2t(ji,jj+1)
1139	!!$ ! ... 2* gradi(v) (v-point i=nibm, time nitm)
1140	!!$ z2dy = ( vsbnd(ji+1,jk,nibm,nitm) - vsbnd(ji-1,jk,nibm,nitm) ) / e1v(ji,jj+1)
1141	!!$ ! ... square of the norm of grad(u)
1142	!!$ z4nor2 = z2dx * z2dx + z2dy * z2dy
1143	!!$ IF( z4nor2 == 0.) THEN
1144	!!$ z4nor2 = 0.000001
1145	!!$ END IF
1146	!!$ ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
1147	!!$ zdt = vsbnd(ji,jk,nibm,nitm2) - vsbnd(ji,jk,nibm,nit)
1148	!!$ ! ... j-phase speed ratio (bounded by -1)
1149	!!$ z05cx = zdt * z2dx / z4nor2
1150	!!$ v_cysbnd(ji,jk)=z05cx*vsmsk(ji,jk)
1151	!!$ END DO
1152	!!$ END DO
1153	!!$ END DO
1154	!!$
1155	!!$ ENDIF
1156	!!$
1157	!!$ END SUBROUTINE obc_rad_south
1158	!!$
1159	!!$#else
1160	!!=================================================================================
1161	!! * MODULE obcrad *
1162	!! Ocean dynamic : Phase velocities for each open boundary
1163	!!=================================================================================
1164	CONTAINS
1165	SUBROUTINE obc_rad( kt ) ! No open boundaries ==> empty routine
1166	INTEGER, INTENT(in) :: kt
1167	WRITE(,) 'obc_rad: You should not have seen this print! error?', kt
1168	END SUBROUTINE obc_rad
1169	#endif
1170
1171	END MODULE obcrad

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

Download in other formats: