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/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/OBC/obcrad.F90 @ 4147

Last change on this file since 4147 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk
  • Property svn:keywords set to Id
File size: 48.1 KB
Line 
1MODULE 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
46CONTAINS
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,jpk*3*3,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,jpk*3*3,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) = tsn(ji  ,jj,jk,jp_tem)*temsk(jj,jk)
218                  tebnd(jj,jk,nibm ,nit) = tsn(ji-1,jj,jk,jp_tem)*temsk(jj,jk)
219                  sebnd(jj,jk,nib  ,nit) = tsn(ji  ,jj,jk,jp_sal)*temsk(jj,jk)
220                  sebnd(jj,jk,nibm ,nit) = tsn(ji-1,jj,jk,jp_sal)*temsk(jj,jk)
221               END DO
222            END DO
223         END DO
224         IF( lk_mpp )   CALL mppobc(tebnd,jpjed,jpjef,jpieob+1,jpk*2*2,2,jpj, numout )
225         IF( lk_mpp )   CALL mppobc(sebnd,jpjed,jpjef,jpieob+1,jpk*2*2,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,jpk*3*3,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,jpk*3*3,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) = tsn(ji   ,jj,jk,jp_tem)*twmsk(jj,jk)
484                  twbnd(jj,jk,nibm ,nit) = tsn(ji+1 ,jj,jk,jp_tem)*twmsk(jj,jk)
485                  swbnd(jj,jk,nib  ,nit) = tsn(ji   ,jj,jk,jp_sal)*twmsk(jj,jk)
486                  swbnd(jj,jk,nibm ,nit) = tsn(ji+1 ,jj,jk,jp_sal)*twmsk(jj,jk)
487               END DO
488            END DO
489         END DO
490         IF( lk_mpp )   CALL mppobc(twbnd,jpjwd,jpjwf,jpiwob,jpk*2*2,2,jpj, numout )
491         IF( lk_mpp )   CALL mppobc(swbnd,jpjwd,jpjwf,jpiwob,jpk*2*2,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,jpk*3*3,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,jpk*3*3,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) = tsn(ji,jj,  jk,jp_tem)*tnmsk(ji,jk)
753                  tnbnd(ji,jk,nibm ,nit) = tsn(ji,jj-1,jk,jp_tem)*tnmsk(ji,jk)
754                  snbnd(ji,jk,nib  ,nit) = tsn(ji,jj,  jk,jp_sal)*tnmsk(ji,jk)
755                  snbnd(ji,jk,nibm ,nit) = tsn(ji,jj-1,jk,jp_sal)*tnmsk(ji,jk)
756               END DO
757            END DO
758         END DO
759         IF( lk_mpp )   CALL mppobc(tnbnd,jpind,jpinf,jpjnob+1,jpk*2*2,1,jpi, numout )
760         IF( lk_mpp )   CALL mppobc(snbnd,jpind,jpinf,jpjnob+1,jpk*2*2,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,jpk*3*3,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,jpk*3*3,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) = tsn(ji,jj   ,jk,jp_tem)*tsmsk(ji,jk)
1025                  tsbnd(ji,jk,nibm ,nit) = tsn(ji,jj+1 ,jk,jp_tem)*tsmsk(ji,jk)
1026                  ssbnd(ji,jk,nib  ,nit) = tsn(ji,jj   ,jk,jp_sal)*tsmsk(ji,jk)
1027                  ssbnd(ji,jk,nibm ,nit) = tsn(ji,jj+1 ,jk,jp_sal)*tsmsk(ji,jk)
1028               END DO
1029            END DO
1030         END DO
1031         IF( lk_mpp )   CALL mppobc(tsbnd,jpisd,jpisf,jpjsob,jpk*2*2,1,jpi, numout )
1032         IF( lk_mpp )   CALL mppobc(ssbnd,jpisd,jpisf,jpjsob,jpk*2*2,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   !!=================================================================================
1164CONTAINS
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
1171END MODULE obcrad
Note: See TracBrowser for help on using the repository browser.