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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC/obcrad.F90 @ 4409

Last change on this file since 4409 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

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