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

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

CT : UPDATE052 : change logical lpXXXobc to lp_obc_XXX for Open Boundaries case

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