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

source: trunk/NEMO/OPA_SRC/OBC/obcrad.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

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