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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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