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.
obcdyn.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 30.3 KB
Line 
1MODULE obcdyn
2#if defined key_obc
3   !!=================================================================================
4   !!                       ***  MODULE  obcdyn  ***
5   !! Ocean dynamics:   Radiation of velocities on each open boundary
6   !!=================================================================================
7
8   !!---------------------------------------------------------------------------------
9   !!   obc_dyn        : call the subroutine for each open boundary
10   !!   obc_dyn_east   : radiation of the east open boundary velocities
11   !!   obc_dyn_west   : radiation of the west open boundary velocities
12   !!   obc_dyn_north  : radiation of the north open boundary velocities
13   !!   obc_dyn_south  : radiation of the south open boundary velocities
14   !!----------------------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE phycst          ! physical constants
21   USE obc_oce         ! ocean open boundary conditions
22   USE lbclnk          ! ???
23   USE lib_mpp         ! ???
24   USE in_out_manager  ! I/O manager
25   USE dynspg_oce
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Accessibility
31   PUBLIC obc_dyn     ! routine called in dynspg_flt (free surface case)
32
33   !! * Module variables
34   INTEGER ::   ji, jj, jk     ! dummy loop indices
35
36   INTEGER ::      & ! ... boundary space indices
37      nib   = 1,   & ! nib   = boundary point
38      nibm  = 2,   & ! nibm  = 1st interior point
39      nibm2 = 3,   & ! nibm2 = 2nd interior point
40                     ! ... boundary time indices
41      nit   = 1,   & ! nit    = now
42      nitm  = 2,   & ! nitm   = before
43      nitm2 = 3      ! nitm2  = before-before
44
45   REAL(wp) ::   rtaue  , rtauw  , rtaun  , rtaus  ,  &
46                 rtauein, rtauwin, rtaunin, rtausin
47
48   !! * Control permutation of array indices
49#  include "oce_ftrans.h90"
50#  include "dom_oce_ftrans.h90"
51#  include "obc_oce_ftrans.h90"
52
53   !!---------------------------------------------------------------------------------
54
55CONTAINS
56
57   SUBROUTINE obc_dyn ( kt )
58      !!------------------------------------------------------------------------------
59      !!                      SUBROUTINE obc_dyn
60      !!                     ********************
61      !! ** Purpose :
62      !!      Compute  dynamics (u,v) at the open boundaries.
63      !!      if defined key_dynspg_flt:
64      !!                 this routine is called by dynspg_flt and updates
65      !!                 ua, va which are the actual velocities (not trends)
66      !!
67      !!      The logical variable lp_obc_east, and/or lp_obc_west, and/or lp_obc_north,
68      !!      and/or lp_obc_south allow the user to determine which boundary is an
69      !!      open one (must be done in the param_obc.h90 file).
70      !!
71      !! ** Reference :
72      !!      Marchesiello P., 1995, these de l'universite J. Fourier, Grenoble, France.
73      !!
74      !! History :
75      !!        !  95-03 (J.-M. Molines) Original, SPEM
76      !!        !  97-07 (G. Madec, J.-M. Molines) addition
77      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90
78      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
79      !!----------------------------------------------------------------------
80      !! * Arguments
81      INTEGER, INTENT( in ) ::   kt
82
83      !!----------------------------------------------------------------------
84      !!  OPA 9.0 , LOCEAN-IPSL (2005)
85      !! $Id$
86      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
87      !!----------------------------------------------------------------------
88
89      ! 0. Local constant initialization
90      ! --------------------------------
91
92      IF( kt == nit000 .OR. ln_rstart) THEN
93         ! ... Boundary restoring coefficient
94         rtaue = 2. * rdt / rdpeob
95         rtauw = 2. * rdt / rdpwob
96         rtaun = 2. * rdt / rdpnob
97         rtaus = 2. * rdt / rdpsob
98         ! ... Boundary restoring coefficient for inflow ( all boundaries)
99         rtauein = 2. * rdt / rdpein 
100         rtauwin = 2. * rdt / rdpwin
101         rtaunin = 2. * rdt / rdpnin
102         rtausin = 2. * rdt / rdpsin 
103      END IF
104
105      IF( lp_obc_east  )   CALL obc_dyn_east ( kt )
106      IF( lp_obc_west  )   CALL obc_dyn_west ( kt )
107      IF( lp_obc_north )   CALL obc_dyn_north( kt )
108      IF( lp_obc_south )   CALL obc_dyn_south( kt )
109
110      IF( lk_mpp ) THEN
111         IF( kt >= nit000+3 .AND. ln_rstart ) THEN
112            CALL lbc_lnk( ub, 'U', -1. )
113            CALL lbc_lnk( vb, 'V', -1. )
114         END IF
115         CALL lbc_lnk( ua, 'U', -1. )
116         CALL lbc_lnk( va, 'V', -1. )
117      ENDIF
118
119   END SUBROUTINE obc_dyn
120
121
122   SUBROUTINE obc_dyn_east ( kt )
123      !!------------------------------------------------------------------------------
124      !!                  ***  SUBROUTINE obc_dyn_east  ***
125      !!             
126      !! ** Purpose :
127      !!      Apply the radiation algorithm on east OBC velocities ua, va using the
128      !!      phase velocities calculated in obc_rad_east subroutine in obcrad.F90 module
129      !!      If the logical lfbceast is .TRUE., there is no radiation but only fixed OBC
130      !!
131      !!  History :
132      !!         ! 95-03 (J.-M. Molines) Original from SPEM
133      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
134      !!         ! 97-12 (M. Imbard) Mpp adaptation
135      !!         ! 00-06 (J.-M. Molines)
136      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
137      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
138      !!------------------------------------------------------------------------------
139      !! * Arguments
140      INTEGER, INTENT( in ) ::   kt
141
142      !! * Local declaration
143      REAL(wp) ::   z05cx, ztau, zin
144      !!------------------------------------------------------------------------------
145
146      ! 1. First three time steps and more if lfbceast is .TRUE.
147      !    In that case open boundary conditions are FIXED.
148      ! --------------------------------------------------------
149
150      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbceast .OR. lk_dynspg_exp ) THEN
151
152         ! 1.1 U zonal velocity   
153         ! --------------------
154#if defined key_z_first
155         DO jj = 1, jpj
156            DO ji = nie0, nie1
157               DO jk = 1, jpkm1
158#else
159         DO ji = nie0, nie1
160            DO jk = 1, jpkm1
161               DO jj = 1, jpj
162#endif
163                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uemsk(jj,jk)) + &
164                                 uemsk(jj,jk)*ufoe(jj,jk)
165               END DO
166            END DO
167         END DO
168
169         ! 1.2 V meridional velocity
170         ! -------------------------
171#if defined key_z_first
172         DO jj = 1, jpj
173            DO ji = nie0+1, nie1+1
174               DO jk = 1, jpkm1
175#else
176         DO ji = nie0+1, nie1+1
177            DO jk = 1, jpkm1
178               DO jj = 1, jpj
179#endif
180                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vemsk(jj,jk)) + &
181                                 vfoe(jj,jk)*vemsk(jj,jk)
182               END DO
183            END DO
184         END DO
185
186      ELSE
187
188      ! 2. Beyond the fourth time step if lfbceast is .FALSE.
189      ! -----------------------------------------------------
190 
191         ! 2.1. u-component of the velocity
192         ! ---------------------------------
193         !
194         !          nibm2      nibm      nib
195         !            |   nibm  |   nib   |///
196         !            |    |    |    |    |///
197         !  jj-line --f----v----f----v----f---
198         !            |    |    |    |    |///
199         !            |         |         |///
200         !  jj-line   u    T    u    T    u///
201         !            |         |         |///
202         !            |    |    |    |    |///
203         !          jpieob-2   jpieob-1   jpieob
204         !                 |         |       
205         !              jpieob-1    jpieob   
206         ! 
207         ! ... If free surface formulation:
208         ! ... radiative conditions on the total part + relaxation toward climatology
209         ! ... (jpjedp1, jpjefm1),jpieob
210#if defined key_z_first
211         DO jj = 1, jpj
212            DO ji = nie0, nie1
213               DO jk = 1, jpkm1
214#else
215         DO ji = nie0, nie1
216            DO jk = 1, jpkm1
217               DO jj = 1, jpj
218#endif
219                  z05cx = u_cxebnd(jj,jk)
220                  z05cx = z05cx / e1t(ji,jj)
221                  z05cx = min( z05cx, 1. )
222         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
223         !           > 0, outflow zin=1, ztau=rtaue
224                  zin = sign( 1., z05cx )
225                  zin = 0.5*( zin + abs(zin) )
226         ! ... for inflow rtauein is used for relaxation coefficient else rtaue
227                  ztau = (1.-zin ) * rtauein  + zin * rtaue
228                  z05cx = z05cx * zin
229         ! ... update ua with radiative or climatological velocity
230                  ua(ji,jj,jk) = ua(ji,jj,jk) * ( 1. - uemsk(jj,jk) ) +          &
231                                 uemsk(jj,jk) * (  ( 1. - z05cx - ztau )         &
232                                 * uebnd(jj,jk,nib ,nitm) + 2.*z05cx               &
233                                 * uebnd(jj,jk,nibm,nit ) + ztau * ufoe (jj,jk) )  &
234                                 / (1. + z05cx)
235               END DO
236            END DO
237         END DO
238
239         ! 2.2 v-component of the velocity
240         ! -------------------------------
241         !
242         !          nibm2       nibm     nib
243         !            |   nibm  |   nib///|///
244         !            |    |    |    |////|///
245         !  jj-line --v----f----v----f----v---
246         !            |    |    |    |////|///
247         !            |    |    |    |////|///
248         !            | jpieob-1 |  jpieob /|///
249         !            |         |         |   
250         !         jpieob-1    jpieob     jpieob+1
251         !
252         ! ... radiative condition
253         ! ... (jpjedp1, jpjefm1), jpieob+1
254#if defined key_z_first
255         DO jj = 1, jpj
256            DO ji = nie0+1, nie1+1
257               DO jk = 1, jpkm1
258#else
259         DO ji = nie0+1, nie1+1
260            DO jk = 1, jpkm1
261               DO jj = 1, jpj
262#endif
263                  z05cx = v_cxebnd(jj,jk) 
264                  z05cx = z05cx / e1f(ji-1,jj)
265                  z05cx = min( z05cx, 1. )
266         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
267         !           > 0, outflow zin=1, ztau=rtaue
268                  zin = sign( 1., z05cx )
269                  zin = 0.5*( zin + abs(zin) )
270         ! ... for inflow rtauein is used for relaxation coefficient else rtaue
271                  ztau = (1.-zin ) * rtauein  + zin * rtaue
272                  z05cx = z05cx * zin
273         ! ... update va with radiative or climatological velocity
274                  va(ji,jj,jk) = va(ji,jj,jk) * (1. - vemsk(jj,jk) ) +          &
275                                 vemsk(jj,jk) * ( ( 1. - z05cx - ztau )         &
276                                 * vebnd(jj,jk,nib ,nitm) + 2.*z05cx              &
277                                 * vebnd(jj,jk,nibm,nit ) + ztau * vfoe(jj,jk) )  &
278                                 / (1. + z05cx)
279               END DO
280            END DO
281         END DO
282
283      END IF
284
285   END SUBROUTINE obc_dyn_east
286
287
288   SUBROUTINE obc_dyn_west ( kt )
289      !!------------------------------------------------------------------------------
290      !!                  ***  SUBROUTINE obc_dyn_west  ***
291      !!                 
292      !! ** Purpose :
293      !!      Apply the radiation algorithm on west OBC velocities ua, va using the
294      !!      phase velocities calculated in obc_rad_west subroutine in obcrad.F90 module
295      !!      If the logical lfbcwest is .TRUE., there is no radiation but only fixed OBC
296      !!
297      !!  History :
298      !!         ! 95-03 (J.-M. Molines) Original from SPEM
299      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
300      !!         ! 97-12 (M. Imbard) Mpp adaptation
301      !!         ! 00-06 (J.-M. Molines)
302      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
303      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
304      !!------------------------------------------------------------------------------
305      !! * Arguments
306      INTEGER, INTENT( in ) ::   kt
307
308      !! * Local declaration
309      REAL(wp) ::   z05cx, ztau, zin
310      !!------------------------------------------------------------------------------
311
312      ! 1. First three time steps and more if lfbcwest is .TRUE.
313      !    In that case open boundary conditions are FIXED.
314      ! --------------------------------------------------------
315
316      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcwest .OR. lk_dynspg_exp ) THEN
317
318         ! 1.1 U zonal velocity
319         ! ---------------------
320#if defined key_z_first
321         DO jj = 1, jpj
322            DO ji = niw0, niw1
323               DO jk = 1, jpkm1
324#else
325         DO ji = niw0, niw1
326            DO jk = 1, jpkm1
327               DO jj = 1, jpj
328#endif
329                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uwmsk(jj,jk)) + &
330                                 uwmsk(jj,jk)*ufow(jj,jk)
331               END DO
332            END DO
333         END DO
334
335         ! 1.2 V meridional velocity
336         ! -------------------------
337#if defined key_z_first
338         DO jj = 1, jpj
339            DO ji = niw0, niw1
340               DO jk = 1, jpkm1
341#else
342         DO ji = niw0, niw1
343            DO jk = 1, jpkm1
344               DO jj = 1, jpj
345#endif
346                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vwmsk(jj,jk)) + &
347                                 vfow(jj,jk)*vwmsk(jj,jk)
348               END DO
349            END DO
350         END DO
351
352      ELSE
353
354      ! 2. Beyond the fourth time step if lfbcwest is .FALSE.
355      ! -----------------------------------------------------
356 
357         ! 2.1. u-component of the velocity
358         ! ---------------------------------
359         !
360         !        nib       nibm     nibm2
361         !      ///|   nib   |   nibm  |
362         !      ///|    |    |    |    |
363         !      ---f----v----f----v----f-- jj-line
364         !      ///|    |    |    |    |
365         !      ///|         |         |
366         !      ///u    T    u    T    u   jj-line
367         !      ///|         |         |
368         !      ///|    |    |    |    |
369         !       jpiwob    jpiwob+1    jpiwob+2
370         !              |         |       
371         !            jpiwob+1    jpiwob+2     
372         !
373         ! ... If free surface formulation:
374         ! ... radiative conditions on the total part + relaxation toward climatology
375         ! ... (jpjwdp1, jpjwfm1), jpiwob
376#if defined key_z_first
377         DO jj = 1, jpj
378            DO ji = niw0, niw1
379               DO jk = 1, jpkm1
380#else
381         DO ji = niw0, niw1
382            DO jk = 1, jpkm1
383               DO jj = 1, jpj
384#endif
385                  z05cx = u_cxwbnd(jj,jk)
386                  z05cx = z05cx / e1t(ji+1,jj)
387                  z05cx = max( z05cx, -1. )
388         ! ... z05c  > 0, inflow  zin=0, ztau=1   
389         !          =< 0, outflow zin=1, ztau=rtauw
390                  zin = sign( 1., -1. * z05cx )
391                  zin = 0.5*( zin + abs(zin) )
392                  ztau = (1.-zin )* rtauwin + zin * rtauw
393                  z05cx = z05cx * zin
394         ! ... update un with radiative or climatological velocity
395                  ua(ji,jj,jk) = ua(ji,jj,jk) * ( 1. - uwmsk(jj,jk) ) +          &
396                                 uwmsk(jj,jk) * ( ( 1. + z05cx - ztau )          &
397                                 * uwbnd(jj,jk,nib ,nitm) - 2.*z05cx               &
398                                 * uwbnd(jj,jk,nibm,nit ) + ztau  * ufow (jj,jk) ) &
399                                 / (1. - z05cx)
400               END DO
401            END DO
402         END DO
403
404         ! 2.2 v-component of the velocity
405         ! -------------------------------
406         !
407         !    nib       nibm     nibm2
408         !  ///|///nib   |   nibm  |  nibm2
409         !  ///|////|    |    |    |    |    |
410         !  ---v----f----v----f----v----f----v-- jj-line
411         !  ///|////|    |    |    |    |    |
412         !  ///|////|    |    |    |    |    |
413         ! jpiwob     jpiwob+1    jpiwob+2
414         !          |         |         |   
415         !        jpiwob   jpiwob+1   jpiwob+2   
416         !
417         ! ... radiative condition plus Raymond-Kuo
418         ! ... (jpjwdp1, jpjwfm1),jpiwob
419#if defined key_z_first
420         DO jj = 1, jpj
421            DO ji = niw0, niw1
422               DO jk = 1, jpkm1
423#else
424         DO ji = niw0, niw1
425            DO jk = 1, jpkm1
426               DO jj = 1, jpj
427#endif
428                  z05cx = v_cxwbnd(jj,jk) 
429                  z05cx = z05cx / e1f(ji,jj)
430                  z05cx = max( z05cx, -1. )
431         ! ... z05cx > 0, inflow  zin=0, ztau=1   
432         !          =< 0, outflow zin=1, ztau=rtauw
433                  zin = sign( 1., -1. * z05cx )
434                  zin = 0.5*( zin + abs(zin) )
435                  ztau = (1.-zin )*rtauwin + zin * rtauw
436                  z05cx = z05cx * zin 
437         ! ... update va with radiative or climatological velocity
438                  va(ji,jj,jk) = va(ji,jj,jk) * (1. - vwmsk(jj,jk) ) +          &
439                                 vwmsk(jj,jk) * ( ( 1. + z05cx - ztau )         &
440                                 * vwbnd(jj,jk,nib ,nitm) - 2.*z05cx              &
441                                 * vwbnd(jj,jk,nibm,nit ) + ztau * vfow (jj,jk) ) &
442                                 / (1. - z05cx)
443                END DO
444             END DO
445         END DO
446
447      END IF
448
449   END SUBROUTINE obc_dyn_west
450
451   SUBROUTINE obc_dyn_north ( kt )
452      !!------------------------------------------------------------------------------
453      !!                     SUBROUTINE obc_dyn_north
454      !!                    *************************
455      !! ** Purpose :
456      !!      Apply the radiation algorithm on north OBC velocities ua, va using the
457      !!      phase velocities calculated in obc_rad_north subroutine in obcrad.F90 module
458      !!      If the logical lfbcnorth is .TRUE., there is no radiation but only fixed OBC
459      !!
460      !!  History :
461      !!         ! 95-03 (J.-M. Molines) Original from SPEM
462      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
463      !!         ! 97-12 (M. Imbard) Mpp adaptation
464      !!         ! 00-06 (J.-M. Molines)
465      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
466      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
467      !!------------------------------------------------------------------------------
468      !! * Arguments
469      INTEGER, INTENT( in ) ::   kt
470
471      !! * Local declaration
472      REAL(wp) ::   z05cx, ztau, zin
473      !!------------------------------------------------------------------------------
474
475      ! 1. First three time steps and more if lfbcnorth is .TRUE.
476      !    In that case open boundary conditions are FIXED.
477      ! ---------------------------------------------------------
478 
479      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcnorth  .OR. lk_dynspg_exp ) THEN
480
481         ! 1.1 U zonal velocity
482         ! --------------------
483         DO jj = njn0+1, njn1+1
484#if defined key_z_first
485            DO ji = 1, jpi
486               DO jk = 1, jpkm1
487#else
488            DO jk = 1, jpkm1
489               DO ji = 1, jpi
490#endif
491                  ua(ji,jj,jk)= ua(ji,jj,jk) * (1.-unmsk(ji,jk)) + &
492                                ufon(ji,jk)*unmsk(ji,jk)
493               END DO
494            END DO
495         END DO
496
497         ! 1.2 V meridional velocity
498         ! -------------------------
499         DO jj = njn0, njn1
500#if defined key_z_first
501            DO ji = 1, jpi
502               DO jk = 1, jpkm1
503#else
504            DO jk = 1, jpkm1
505               DO ji = 1, jpi
506#endif
507                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vnmsk(ji,jk)) + &
508                                vfon(ji,jk)*vnmsk(ji,jk)
509               END DO
510            END DO
511         END DO
512
513      ELSE
514
515      ! 2. Beyond the fourth time step if lfbcnorth is .FALSE.
516      ! ------------------------------------------------------
517
518         ! 2.1. u-component of the velocity
519         ! --------------------------------
520         !
521         !            ji-row
522         !              |
523         !       nib ///u//////  jpjnob + 1
524         !         /////|//////
525         !     nib -----f-----   jpjnob
526         !              |   
527         !      nibm--  u   ---- jpjnob
528         !              |       
529         !    nibm -----f-----   jpjnob-1
530         !              |       
531         !     nibm2--  u   ---- jpjnob-1
532         !              |       
533         !   nibm2 -----f-----   jpjnob-2
534         !              |
535         !
536         ! ... radiative condition
537         ! ... jpjnob+1,(jpindp1, jpinfm1)
538         DO jj = njn0+1, njn1+1
539#if defined key_z_first
540            DO ji = 1, jpi
541               DO jk = 1, jpkm1
542#else
543            DO jk = 1, jpkm1
544               DO ji = 1, jpi
545#endif
546                  z05cx= u_cynbnd(ji,jk) 
547                  z05cx = z05cx / e2f(ji, jj-1)
548                  z05cx = min( z05cx, 1. )
549         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
550         !           > 0, outflow zin=1, ztau=rtaun
551                  zin = sign( 1., z05cx )
552                  zin = 0.5*( zin + abs(zin) )
553         ! ... for inflow rtaunin is used for relaxation coefficient else rtaun
554                  ztau = (1.-zin ) * rtaunin  + zin * rtaun
555         ! ... for u, when inflow, ufon is prescribed
556                  z05cx = z05cx * zin
557         ! ... update un with radiative or climatological velocity
558                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-unmsk(ji,jk)) +             &
559                                 unmsk(ji,jk) * ( ( 1. - z05cx - ztau )         &
560                                 * unbnd(ji,jk,nib ,nitm) + 2.*z05cx              &
561                                 * unbnd(ji,jk,nibm,nit ) + ztau * ufon (ji,jk) ) &
562                                 / (1. + z05cx)
563               END DO
564            END DO
565         END DO
566
567         ! 2.2 v-component of the velocity
568         ! -------------------------------
569         !
570         !                ji-row    ji-row
571         !              |         |
572         !         /////|/////////////////
573         !    nib  -----f----v----f----  jpjnob
574         !              |         |
575         !      nib  -  u -- T -- u ---- jpjnob
576         !              |         |
577         !   nibm  -----f----v----f----  jpjnob-1
578         !              |         |
579         !     nibm --  u -- T -- u ---  jpjnob-1
580         !              |         |
581         !   nibm2 -----f----v----f----  jpjnob-2
582         !              |         |
583         !
584         ! ... Free surface formulation:
585         ! ... radiative conditions on the total part + relaxation toward climatology
586         ! ... jpjnob,(jpindp1, jpinfm1)
587         DO jj = njn0, njn1
588#if defined key_z_first
589            DO ji = 1, jpi
590               DO jk = 1, jpkm1
591#else
592            DO jk = 1, jpkm1
593               DO ji = 1, jpi
594#endif
595         ! ... 2* gradj(v) (T-point i=nibm, time mean)
596                  z05cx = v_cynbnd(ji,jk)
597                  z05cx = z05cx / e2t(ji,jj)
598                  z05cx = min( z05cx, 1. )
599         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
600         !           > 0, outflow zin=1, ztau=rtaun
601                  zin = sign( 1., z05cx )
602                  zin = 0.5*( zin + abs(zin) )
603         ! ... for inflow rtaunin is used for relaxation coefficient else rtaun
604                  ztau = (1.-zin ) * rtaunin + zin * rtaun
605                  z05cx = z05cx * zin
606         ! ... update va with radiative or climatological velocity
607                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vnmsk(ji,jk)) +             &
608                                 vnmsk(ji,jk) * ( ( 1. - z05cx - ztau )         &
609                                 * vnbnd(ji,jk,nib ,nitm) + 2.*z05cx              &
610                                 * vnbnd(ji,jk,nibm,nit ) + ztau * vfon (ji,jk) ) &
611                                 / (1. + z05cx)
612               END DO
613            END DO
614         END DO
615      END IF
616
617   END SUBROUTINE obc_dyn_north
618
619   SUBROUTINE obc_dyn_south ( kt )
620      !!------------------------------------------------------------------------------
621      !!                     SUBROUTINE obc_dyn_south
622      !!                    *************************
623      !! ** Purpose :
624      !!      Apply the radiation algorithm on south OBC velocities ua, va using the
625      !!      phase velocities calculated in obc_rad_south subroutine in obcrad.F90 module
626      !!      If the logical lfbcsouth is .TRUE., there is no radiation but only fixed OBC
627      !!
628      !!  History :
629      !!         ! 95-03 (J.-M. Molines) Original from SPEM
630      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
631      !!         ! 97-12 (M. Imbard) Mpp adaptation
632      !!         ! 00-06 (J.-M. Molines)
633      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
634      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
635      !!------------------------------------------------------------------------------
636      !! * Arguments
637      INTEGER, INTENT( in ) ::   kt
638
639      !! * Local declaration
640      REAL(wp) ::   z05cx, ztau, zin
641
642      !!------------------------------------------------------------------------------
643      !!  OPA 8.5, LODYC-IPSL (2002)
644      !!------------------------------------------------------------------------------
645
646      ! 1. First three time steps and more if lfbcsouth is .TRUE.
647      !    In that case open boundary conditions are FIXED.
648      ! ---------------------------------------------------------
649
650      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcsouth  .OR. lk_dynspg_exp ) THEN
651
652         ! 1.1 U zonal velocity
653         ! --------------------
654         DO jj = njs0, njs1
655#if defined key_z_first
656            DO ji = 1, jpi
657               DO jk = 1, jpkm1
658#else
659            DO jk = 1, jpkm1
660               DO ji = 1, jpi
661#endif
662                  ua(ji,jj,jk)= ua(ji,jj,jk) * (1.-usmsk(ji,jk)) + &
663                                usmsk(ji,jk) * ufos(ji,jk)
664               END DO
665            END DO
666         END DO
667
668         ! 1.2 V meridional velocity
669         ! -------------------------
670         DO jj = njs0, njs1
671#if defined key_z_first
672            DO ji = 1, jpi
673               DO jk = 1, jpkm1
674#else
675            DO jk = 1, jpkm1
676               DO ji = 1, jpi
677#endif
678                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vsmsk(ji,jk)) + &
679                                vsmsk(ji,jk) * vfos(ji,jk)
680               END DO
681            END DO
682         END DO
683
684      ELSE
685
686      ! 2. Beyond the fourth time step if lfbcsouth is .FALSE.
687      ! ------------------------------------------------------
688
689         ! 2.1. u-component of the velocity
690         ! --------------------------------
691         !
692         !            ji-row
693         !              |
694         !   nibm2 -----f-----   jpjsob +2
695         !              |   
696         !    nibm2 --  u   ---- jpjsob +2
697         !              |       
698         !    nibm -----f-----   jpjsob +1
699         !              |       
700         !    nibm  --  u   ---- jpjsob +1
701         !              |       
702         !    nib  -----f-----   jpjsob
703         !         /////|//////
704         !    nib   ////u/////   jpjsob
705         !
706         ! ... radiative condition plus Raymond-Kuo
707         ! ... jpjsob,(jpisdp1, jpisfm1)
708         DO jj = njs0, njs1
709#if defined key_z_first
710            DO ji = 1, jpi
711               DO jk = 1, jpkm1
712#else
713            DO jk = 1, jpkm1
714               DO ji = 1, jpi
715#endif
716                  z05cx= u_cysbnd(ji,jk) 
717                  z05cx = z05cx / e2f(ji, jj)
718                  z05cx = max( z05cx, -1. )
719         ! ... z05cx > 0, inflow  zin=0, ztau=1
720         !          =< 0, outflow zin=1, ztau=rtaus
721                  zin = sign( 1., -1. * z05cx )
722                  zin = 0.5*( zin + abs(zin) )
723                  ztau = (1.-zin ) * rtausin + zin * rtaus
724                  z05cx = z05cx * zin 
725         ! ... update ua with radiative or climatological velocity
726                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-usmsk(ji,jk)) +              &
727                                 usmsk(ji,jk) * (  ( 1. + z05cx - ztau )         &
728                                 * usbnd(ji,jk,nib ,nitm) - 2.*z05cx               &
729                                 * usbnd(ji,jk,nibm,nit ) + ztau * ufos (ji,jk) )  &
730                                 / (1. - z05cx)
731               END DO
732            END DO
733         END DO
734
735         ! 2.2 v-component of the velocity
736         ! -------------------------------
737         !
738         !                ji-row    ji-row
739         !              |         |
740         !  nibm2  -----f----v----f----  jpjsob+2
741         !              |         |
742         !    nibm   -  u -- T -- u ---- jpjsob+2
743         !              |         |
744         !   nibm  -----f----v----f----  jpjsob+1
745         !              |         |
746         !   nib    --  u -- T -- u ---  jpjsob+1
747         !              |         |
748         !   nib   -----f----v----f----  jpjsob
749         !         /////////////////////
750         !
751         ! ... Free surface formulation:
752         ! ... radiative conditions on the total part + relaxation toward climatology
753         ! ... jpjsob,(jpisdp1,jpisfm1)
754         DO jj = njs0, njs1
755#if defined key_z_first
756            DO ji = 1, jpi
757               DO jk = 1, jpkm1
758#else
759            DO jk = 1, jpkm1
760               DO ji = 1, jpi
761#endif
762                  z05cx = v_cysbnd(ji,jk)
763                  z05cx = z05cx / e2t(ji,jj+1)
764                  z05cx = max( z05cx, -1. )
765         ! ... z05c > 0, inflow  zin=0, ztau=1
766         !         =< 0, outflow zin=1, ztau=rtaus
767                  zin = sign( 1., -1. * z05cx )
768                  zin = 0.5*( zin + abs(zin) )
769                  ztau = (1.-zin )*rtausin + zin * rtaus
770                  z05cx = z05cx * zin
771         ! ... update va with radiative or climatological velocity
772                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vsmsk(ji,jk)) +             &
773                                 vsmsk(ji,jk) * ( ( 1. + z05cx - ztau )         &
774                                 * vsbnd(ji,jk,nib ,nitm) - 2.*z05cx              &
775                                 * vsbnd(ji,jk,nibm,nit ) + ztau * vfos (ji,jk) ) &
776                                 / (1. - z05cx)
777               END DO
778            END DO
779         END DO
780      END IF
781
782   END SUBROUTINE obc_dyn_south
783#else
784   !!=================================================================================
785   !!                       ***  MODULE  obcdyn  ***
786   !! Ocean dynamics:   Radiation of velocities on each open boundary
787   !!=================================================================================
788CONTAINS
789
790   SUBROUTINE obc_dyn
791                              ! No open boundaries ==> empty routine
792   END SUBROUTINE obc_dyn
793#endif
794
795END MODULE obcdyn
Note: See TracBrowser for help on using the repository browser.