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

source: trunk/NEMO/OPA_SRC/OBC/obcdyn.F90 @ 1152

Last change on this file since 1152 was 1152, checked in by rblod, 16 years ago

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