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

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

CT : UPDATE001 : First major NEMO update

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