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

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBC/obcrad.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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