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 @ 78

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

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