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

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

Initial revision

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