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 branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90 @ 2888

Last change on this file since 2888 was 2888, checked in by davestorkey, 13 years ago

Move changes into updated BDY module and restore old OBC code.
(Full merge to take place next year).
File size: 28.5 KB
Line 
1MODULE obcdyn
2#if defined key_obc
3   !!=================================================================================
4   !!                       ***  MODULE  obcdyn  ***
5   !! Ocean dynamics:   Radiation of velocities on each open boundary
6   !!=================================================================================
7
8   !!---------------------------------------------------------------------------------
9   !!   obc_dyn        : call the subroutine for each open boundary
10   !!   obc_dyn_east   : radiation of the east open boundary velocities
11   !!   obc_dyn_west   : radiation of the west open boundary velocities
12   !!   obc_dyn_north  : radiation of the north open boundary velocities
13   !!   obc_dyn_south  : radiation of the south open boundary velocities
14   !!----------------------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------------------
17   !! * Modules used
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE phycst          ! physical constants
21   USE obc_oce         ! ocean open boundary conditions
22   USE lbclnk          ! ???
23   USE lib_mpp         ! ???
24   USE in_out_manager  ! I/O manager
25   USE dynspg_oce
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Accessibility
31   PUBLIC obc_dyn     ! routine called in dynspg_flt (free surface 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_flt:
59      !!                 this routine is called by dynspg_flt and updates
60      !!                 ua, va which are the actual velocities (not trends)
61      !!
62      !!      The logical variable lp_obc_east, and/or lp_obc_west, and/or lp_obc_north,
63      !!      and/or lp_obc_south allow the user to determine which boundary is an
64      !!      open one (must be done in the param_obc.h90 file).
65      !!
66      !! ** Reference :
67      !!      Marchesiello P., 1995, these de l'universite J. Fourier, Grenoble, France.
68      !!
69      !! History :
70      !!        !  95-03 (J.-M. Molines) Original, SPEM
71      !!        !  97-07 (G. Madec, J.-M. Molines) addition
72      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90
73      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
74      !!----------------------------------------------------------------------
75      !! * Arguments
76      INTEGER, INTENT( in ) ::   kt
77
78      !!----------------------------------------------------------------------
79      !!  OPA 9.0 , LOCEAN-IPSL (2005)
80      !! $Id: obcdyn.F90 1528 2009-07-23 14:38:47Z rblod $
81      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
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      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
133      !!------------------------------------------------------------------------------
134      !! * Arguments
135      INTEGER, INTENT( in ) ::   kt
136
137      !! * Local declaration
138      REAL(wp) ::   z05cx, ztau, zin
139      !!------------------------------------------------------------------------------
140
141      ! 1. First three time steps and more if lfbceast is .TRUE.
142      !    In that case open boundary conditions are FIXED.
143      ! --------------------------------------------------------
144
145      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbceast .OR. lk_dynspg_exp ) THEN
146
147         ! 1.1 U zonal velocity   
148         ! --------------------
149         DO ji = nie0, nie1
150            DO jk = 1, jpkm1
151               DO jj = 1, jpj
152                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uemsk(jj,jk)) + &
153                                 uemsk(jj,jk)*ufoe(jj,jk)
154               END DO
155            END DO
156         END DO
157
158         ! 1.2 V meridional velocity
159         ! -------------------------
160         DO ji = nie0+1, nie1+1
161            DO jk = 1, jpkm1
162               DO jj = 1, jpj
163                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vemsk(jj,jk)) + &
164                                 vfoe(jj,jk)*vemsk(jj,jk)
165               END DO
166            END DO
167         END DO
168
169      ELSE
170
171      ! 2. Beyond the fourth time step if lfbceast is .FALSE.
172      ! -----------------------------------------------------
173 
174         ! 2.1. u-component of the velocity
175         ! ---------------------------------
176         !
177         !          nibm2      nibm      nib
178         !            |   nibm  |   nib   |///
179         !            |    |    |    |    |///
180         !  jj-line --f----v----f----v----f---
181         !            |    |    |    |    |///
182         !            |         |         |///
183         !  jj-line   u    T    u    T    u///
184         !            |         |         |///
185         !            |    |    |    |    |///
186         !          jpieob-2   jpieob-1   jpieob
187         !                 |         |       
188         !              jpieob-1    jpieob   
189         ! 
190         ! ... If free surface formulation:
191         ! ... radiative conditions on the total part + relaxation toward climatology
192         ! ... (jpjedp1, jpjefm1),jpieob
193         DO ji = nie0, nie1
194            DO jk = 1, jpkm1
195               DO jj = 1, jpj
196                  z05cx = u_cxebnd(jj,jk)
197                  z05cx = z05cx / e1t(ji,jj)
198                  z05cx = min( z05cx, 1. )
199         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
200         !           > 0, outflow zin=1, ztau=rtaue
201                  zin = sign( 1., z05cx )
202                  zin = 0.5*( zin + abs(zin) )
203         ! ... for inflow rtauein is used for relaxation coefficient else rtaue
204                  ztau = (1.-zin ) * rtauein  + zin * rtaue
205                  z05cx = z05cx * zin
206         ! ... update ua with radiative or climatological velocity
207                  ua(ji,jj,jk) = ua(ji,jj,jk) * ( 1. - uemsk(jj,jk) ) +          &
208                                 uemsk(jj,jk) * (  ( 1. - z05cx - ztau )         &
209                                 * uebnd(jj,jk,nib ,nitm) + 2.*z05cx               &
210                                 * uebnd(jj,jk,nibm,nit ) + ztau * ufoe (jj,jk) )  &
211                                 / (1. + z05cx)
212               END DO
213            END DO
214         END DO
215
216         ! 2.2 v-component of the velocity
217         ! -------------------------------
218         !
219         !          nibm2       nibm     nib
220         !            |   nibm  |   nib///|///
221         !            |    |    |    |////|///
222         !  jj-line --v----f----v----f----v---
223         !            |    |    |    |////|///
224         !            |    |    |    |////|///
225         !            | jpieob-1 |  jpieob /|///
226         !            |         |         |   
227         !         jpieob-1    jpieob     jpieob+1
228         !
229         ! ... radiative condition
230         ! ... (jpjedp1, jpjefm1), jpieob+1
231         DO ji = nie0+1, nie1+1
232            DO jk = 1, jpkm1
233               DO jj = 1, jpj
234                  z05cx = v_cxebnd(jj,jk) 
235                  z05cx = z05cx / e1f(ji-1,jj)
236                  z05cx = min( z05cx, 1. )
237         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
238         !           > 0, outflow zin=1, ztau=rtaue
239                  zin = sign( 1., z05cx )
240                  zin = 0.5*( zin + abs(zin) )
241         ! ... for inflow rtauein is used for relaxation coefficient else rtaue
242                  ztau = (1.-zin ) * rtauein  + zin * rtaue
243                  z05cx = z05cx * zin
244         ! ... update va with radiative or climatological velocity
245                  va(ji,jj,jk) = va(ji,jj,jk) * (1. - vemsk(jj,jk) ) +          &
246                                 vemsk(jj,jk) * ( ( 1. - z05cx - ztau )         &
247                                 * vebnd(jj,jk,nib ,nitm) + 2.*z05cx              &
248                                 * vebnd(jj,jk,nibm,nit ) + ztau * vfoe(jj,jk) )  &
249                                 / (1. + z05cx)
250               END DO
251            END DO
252         END DO
253
254      END IF
255
256   END SUBROUTINE obc_dyn_east
257
258
259   SUBROUTINE obc_dyn_west ( kt )
260      !!------------------------------------------------------------------------------
261      !!                  ***  SUBROUTINE obc_dyn_west  ***
262      !!                 
263      !! ** Purpose :
264      !!      Apply the radiation algorithm on west OBC velocities ua, va using the
265      !!      phase velocities calculated in obc_rad_west subroutine in obcrad.F90 module
266      !!      If the logical lfbcwest is .TRUE., there is no radiation but only fixed OBC
267      !!
268      !!  History :
269      !!         ! 95-03 (J.-M. Molines) Original from SPEM
270      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
271      !!         ! 97-12 (M. Imbard) Mpp adaptation
272      !!         ! 00-06 (J.-M. Molines)
273      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
274      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
275      !!------------------------------------------------------------------------------
276      !! * Arguments
277      INTEGER, INTENT( in ) ::   kt
278
279      !! * Local declaration
280      REAL(wp) ::   z05cx, ztau, zin
281      !!------------------------------------------------------------------------------
282
283      ! 1. First three time steps and more if lfbcwest is .TRUE.
284      !    In that case open boundary conditions are FIXED.
285      ! --------------------------------------------------------
286
287      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcwest .OR. lk_dynspg_exp ) THEN
288
289         ! 1.1 U zonal velocity
290         ! ---------------------
291         DO ji = niw0, niw1
292            DO jk = 1, jpkm1
293               DO jj = 1, jpj
294                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-uwmsk(jj,jk)) + &
295                                 uwmsk(jj,jk)*ufow(jj,jk)
296               END DO
297            END DO
298         END DO
299
300         ! 1.2 V meridional velocity
301         ! -------------------------
302         DO ji = niw0, niw1
303            DO jk = 1, jpkm1
304               DO jj = 1, jpj
305                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vwmsk(jj,jk)) + &
306                                 vfow(jj,jk)*vwmsk(jj,jk)
307               END DO
308            END DO
309         END DO
310
311      ELSE
312
313      ! 2. Beyond the fourth time step if lfbcwest is .FALSE.
314      ! -----------------------------------------------------
315 
316         ! 2.1. u-component of the velocity
317         ! ---------------------------------
318         !
319         !        nib       nibm     nibm2
320         !      ///|   nib   |   nibm  |
321         !      ///|    |    |    |    |
322         !      ---f----v----f----v----f-- jj-line
323         !      ///|    |    |    |    |
324         !      ///|         |         |
325         !      ///u    T    u    T    u   jj-line
326         !      ///|         |         |
327         !      ///|    |    |    |    |
328         !       jpiwob    jpiwob+1    jpiwob+2
329         !              |         |       
330         !            jpiwob+1    jpiwob+2     
331         !
332         ! ... If free surface formulation:
333         ! ... radiative conditions on the total part + relaxation toward climatology
334         ! ... (jpjwdp1, jpjwfm1), jpiwob
335         DO ji = niw0, niw1
336            DO jk = 1, jpkm1
337               DO jj = 1, jpj
338                  z05cx = u_cxwbnd(jj,jk)
339                  z05cx = z05cx / e1t(ji+1,jj)
340                  z05cx = max( z05cx, -1. )
341         ! ... z05c  > 0, inflow  zin=0, ztau=1   
342         !          =< 0, outflow zin=1, ztau=rtauw
343                  zin = sign( 1., -1. * z05cx )
344                  zin = 0.5*( zin + abs(zin) )
345                  ztau = (1.-zin )* rtauwin + zin * rtauw
346                  z05cx = z05cx * zin
347         ! ... update un with radiative or climatological velocity
348                  ua(ji,jj,jk) = ua(ji,jj,jk) * ( 1. - uwmsk(jj,jk) ) +          &
349                                 uwmsk(jj,jk) * ( ( 1. + z05cx - ztau )          &
350                                 * uwbnd(jj,jk,nib ,nitm) - 2.*z05cx               &
351                                 * uwbnd(jj,jk,nibm,nit ) + ztau  * ufow (jj,jk) ) &
352                                 / (1. - z05cx)
353               END DO
354            END DO
355         END DO
356
357         ! 2.2 v-component of the velocity
358         ! -------------------------------
359         !
360         !    nib       nibm     nibm2
361         !  ///|///nib   |   nibm  |  nibm2
362         !  ///|////|    |    |    |    |    |
363         !  ---v----f----v----f----v----f----v-- jj-line
364         !  ///|////|    |    |    |    |    |
365         !  ///|////|    |    |    |    |    |
366         ! jpiwob     jpiwob+1    jpiwob+2
367         !          |         |         |   
368         !        jpiwob   jpiwob+1   jpiwob+2   
369         !
370         ! ... radiative condition plus Raymond-Kuo
371         ! ... (jpjwdp1, jpjwfm1),jpiwob
372         DO ji = niw0, niw1
373            DO jk = 1, jpkm1
374               DO jj = 1, jpj
375                  z05cx = v_cxwbnd(jj,jk) 
376                  z05cx = z05cx / e1f(ji,jj)
377                  z05cx = max( z05cx, -1. )
378         ! ... z05cx > 0, inflow  zin=0, ztau=1   
379         !          =< 0, outflow zin=1, ztau=rtauw
380                  zin = sign( 1., -1. * z05cx )
381                  zin = 0.5*( zin + abs(zin) )
382                  ztau = (1.-zin )*rtauwin + zin * rtauw
383                  z05cx = z05cx * zin 
384         ! ... update va with radiative or climatological velocity
385                  va(ji,jj,jk) = va(ji,jj,jk) * (1. - vwmsk(jj,jk) ) +          &
386                                 vwmsk(jj,jk) * ( ( 1. + z05cx - ztau )         &
387                                 * vwbnd(jj,jk,nib ,nitm) - 2.*z05cx              &
388                                 * vwbnd(jj,jk,nibm,nit ) + ztau * vfow (jj,jk) ) &
389                                 / (1. - z05cx)
390                END DO
391             END DO
392         END DO
393
394      END IF
395
396   END SUBROUTINE obc_dyn_west
397
398   SUBROUTINE obc_dyn_north ( kt )
399      !!------------------------------------------------------------------------------
400      !!                     SUBROUTINE obc_dyn_north
401      !!                    *************************
402      !! ** Purpose :
403      !!      Apply the radiation algorithm on north OBC velocities ua, va using the
404      !!      phase velocities calculated in obc_rad_north subroutine in obcrad.F90 module
405      !!      If the logical lfbcnorth is .TRUE., there is no radiation but only fixed OBC
406      !!
407      !!  History :
408      !!         ! 95-03 (J.-M. Molines) Original from SPEM
409      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
410      !!         ! 97-12 (M. Imbard) Mpp adaptation
411      !!         ! 00-06 (J.-M. Molines)
412      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
413      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
414      !!------------------------------------------------------------------------------
415      !! * Arguments
416      INTEGER, INTENT( in ) ::   kt
417
418      !! * Local declaration
419      REAL(wp) ::   z05cx, ztau, zin
420      !!------------------------------------------------------------------------------
421
422      ! 1. First three time steps and more if lfbcnorth is .TRUE.
423      !    In that case open boundary conditions are FIXED.
424      ! ---------------------------------------------------------
425 
426      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcnorth  .OR. lk_dynspg_exp ) THEN
427
428         ! 1.1 U zonal velocity
429         ! --------------------
430         DO jj = njn0+1, njn1+1
431            DO jk = 1, jpkm1
432               DO ji = 1, jpi
433                  ua(ji,jj,jk)= ua(ji,jj,jk) * (1.-unmsk(ji,jk)) + &
434                                ufon(ji,jk)*unmsk(ji,jk)
435               END DO
436            END DO
437         END DO
438
439         ! 1.2 V meridional velocity
440         ! -------------------------
441         DO jj = njn0, njn1
442            DO jk = 1, jpkm1
443               DO ji = 1, jpi
444                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vnmsk(ji,jk)) + &
445                                vfon(ji,jk)*vnmsk(ji,jk)
446               END DO
447            END DO
448         END DO
449
450      ELSE
451
452      ! 2. Beyond the fourth time step if lfbcnorth is .FALSE.
453      ! ------------------------------------------------------
454
455         ! 2.1. u-component of the velocity
456         ! --------------------------------
457         !
458         !            ji-row
459         !              |
460         !       nib ///u//////  jpjnob + 1
461         !         /////|//////
462         !     nib -----f-----   jpjnob
463         !              |   
464         !      nibm--  u   ---- jpjnob
465         !              |       
466         !    nibm -----f-----   jpjnob-1
467         !              |       
468         !     nibm2--  u   ---- jpjnob-1
469         !              |       
470         !   nibm2 -----f-----   jpjnob-2
471         !              |
472         !
473         ! ... radiative condition
474         ! ... jpjnob+1,(jpindp1, jpinfm1)
475         DO jj = njn0+1, njn1+1
476            DO jk = 1, jpkm1
477               DO ji = 1, jpi
478                  z05cx= u_cynbnd(ji,jk) 
479                  z05cx = z05cx / e2f(ji, jj-1)
480                  z05cx = min( z05cx, 1. )
481         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
482         !           > 0, outflow zin=1, ztau=rtaun
483                  zin = sign( 1., z05cx )
484                  zin = 0.5*( zin + abs(zin) )
485         ! ... for inflow rtaunin is used for relaxation coefficient else rtaun
486                  ztau = (1.-zin ) * rtaunin  + zin * rtaun
487         ! ... for u, when inflow, ufon is prescribed
488                  z05cx = z05cx * zin
489         ! ... update un with radiative or climatological velocity
490                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-unmsk(ji,jk)) +             &
491                                 unmsk(ji,jk) * ( ( 1. - z05cx - ztau )         &
492                                 * unbnd(ji,jk,nib ,nitm) + 2.*z05cx              &
493                                 * unbnd(ji,jk,nibm,nit ) + ztau * ufon (ji,jk) ) &
494                                 / (1. + z05cx)
495               END DO
496            END DO
497         END DO
498
499         ! 2.2 v-component of the velocity
500         ! -------------------------------
501         !
502         !                ji-row    ji-row
503         !              |         |
504         !         /////|/////////////////
505         !    nib  -----f----v----f----  jpjnob
506         !              |         |
507         !      nib  -  u -- T -- u ---- jpjnob
508         !              |         |
509         !   nibm  -----f----v----f----  jpjnob-1
510         !              |         |
511         !     nibm --  u -- T -- u ---  jpjnob-1
512         !              |         |
513         !   nibm2 -----f----v----f----  jpjnob-2
514         !              |         |
515         !
516         ! ... Free surface formulation:
517         ! ... radiative conditions on the total part + relaxation toward climatology
518         ! ... jpjnob,(jpindp1, jpinfm1)
519         DO jj = njn0, njn1
520            DO jk = 1, jpkm1
521               DO ji = 1, jpi
522         ! ... 2* gradj(v) (T-point i=nibm, time mean)
523                  z05cx = v_cynbnd(ji,jk)
524                  z05cx = z05cx / e2t(ji,jj)
525                  z05cx = min( z05cx, 1. )
526         ! ... z05cx=< 0, inflow  zin=0, ztau=1   
527         !           > 0, outflow zin=1, ztau=rtaun
528                  zin = sign( 1., z05cx )
529                  zin = 0.5*( zin + abs(zin) )
530         ! ... for inflow rtaunin is used for relaxation coefficient else rtaun
531                  ztau = (1.-zin ) * rtaunin + zin * rtaun
532                  z05cx = z05cx * zin
533         ! ... update va with radiative or climatological velocity
534                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vnmsk(ji,jk)) +             &
535                                 vnmsk(ji,jk) * ( ( 1. - z05cx - ztau )         &
536                                 * vnbnd(ji,jk,nib ,nitm) + 2.*z05cx              &
537                                 * vnbnd(ji,jk,nibm,nit ) + ztau * vfon (ji,jk) ) &
538                                 / (1. + z05cx)
539               END DO
540            END DO
541         END DO
542      END IF
543
544   END SUBROUTINE obc_dyn_north
545
546   SUBROUTINE obc_dyn_south ( kt )
547      !!------------------------------------------------------------------------------
548      !!                     SUBROUTINE obc_dyn_south
549      !!                    *************************
550      !! ** Purpose :
551      !!      Apply the radiation algorithm on south OBC velocities ua, va using the
552      !!      phase velocities calculated in obc_rad_south subroutine in obcrad.F90 module
553      !!      If the logical lfbcsouth is .TRUE., there is no radiation but only fixed OBC
554      !!
555      !!  History :
556      !!         ! 95-03 (J.-M. Molines) Original from SPEM
557      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
558      !!         ! 97-12 (M. Imbard) Mpp adaptation
559      !!         ! 00-06 (J.-M. Molines)
560      !!    8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
561      !!    9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
562      !!------------------------------------------------------------------------------
563      !! * Arguments
564      INTEGER, INTENT( in ) ::   kt
565
566      !! * Local declaration
567      REAL(wp) ::   z05cx, ztau, zin
568
569      !!------------------------------------------------------------------------------
570      !!  OPA 8.5, LODYC-IPSL (2002)
571      !!------------------------------------------------------------------------------
572
573      ! 1. First three time steps and more if lfbcsouth is .TRUE.
574      !    In that case open boundary conditions are FIXED.
575      ! ---------------------------------------------------------
576
577      IF( ( kt < nit000+3 .AND. .NOT.ln_rstart ) .OR. lfbcsouth  .OR. lk_dynspg_exp ) THEN
578
579         ! 1.1 U zonal velocity
580         ! --------------------
581         DO jj = njs0, njs1
582            DO jk = 1, jpkm1
583               DO ji = 1, jpi
584                  ua(ji,jj,jk)= ua(ji,jj,jk) * (1.-usmsk(ji,jk)) + &
585                                usmsk(ji,jk) * ufos(ji,jk)
586               END DO
587            END DO
588         END DO
589
590         ! 1.2 V meridional velocity
591         ! -------------------------
592         DO jj = njs0, njs1
593            DO jk = 1, jpkm1
594               DO ji = 1, jpi
595                  va(ji,jj,jk)= va(ji,jj,jk) * (1.-vsmsk(ji,jk)) + &
596                                vsmsk(ji,jk) * vfos(ji,jk)
597               END DO
598            END DO
599         END DO
600
601      ELSE
602
603      ! 2. Beyond the fourth time step if lfbcsouth is .FALSE.
604      ! ------------------------------------------------------
605
606         ! 2.1. u-component of the velocity
607         ! --------------------------------
608         !
609         !            ji-row
610         !              |
611         !   nibm2 -----f-----   jpjsob +2
612         !              |   
613         !    nibm2 --  u   ---- jpjsob +2
614         !              |       
615         !    nibm -----f-----   jpjsob +1
616         !              |       
617         !    nibm  --  u   ---- jpjsob +1
618         !              |       
619         !    nib  -----f-----   jpjsob
620         !         /////|//////
621         !    nib   ////u/////   jpjsob
622         !
623         ! ... radiative condition plus Raymond-Kuo
624         ! ... jpjsob,(jpisdp1, jpisfm1)
625         DO jj = njs0, njs1
626            DO jk = 1, jpkm1
627               DO ji = 1, jpi
628                  z05cx= u_cysbnd(ji,jk) 
629                  z05cx = z05cx / e2f(ji, jj)
630                  z05cx = max( z05cx, -1. )
631         ! ... z05cx > 0, inflow  zin=0, ztau=1
632         !          =< 0, outflow zin=1, ztau=rtaus
633                  zin = sign( 1., -1. * z05cx )
634                  zin = 0.5*( zin + abs(zin) )
635                  ztau = (1.-zin ) * rtausin + zin * rtaus
636                  z05cx = z05cx * zin 
637         ! ... update ua with radiative or climatological velocity
638                  ua(ji,jj,jk) = ua(ji,jj,jk) * (1.-usmsk(ji,jk)) +              &
639                                 usmsk(ji,jk) * (  ( 1. + z05cx - ztau )         &
640                                 * usbnd(ji,jk,nib ,nitm) - 2.*z05cx               &
641                                 * usbnd(ji,jk,nibm,nit ) + ztau * ufos (ji,jk) )  &
642                                 / (1. - z05cx)
643               END DO
644            END DO
645         END DO
646
647         ! 2.2 v-component of the velocity
648         ! -------------------------------
649         !
650         !                ji-row    ji-row
651         !              |         |
652         !  nibm2  -----f----v----f----  jpjsob+2
653         !              |         |
654         !    nibm   -  u -- T -- u ---- jpjsob+2
655         !              |         |
656         !   nibm  -----f----v----f----  jpjsob+1
657         !              |         |
658         !   nib    --  u -- T -- u ---  jpjsob+1
659         !              |         |
660         !   nib   -----f----v----f----  jpjsob
661         !         /////////////////////
662         !
663         ! ... Free surface formulation:
664         ! ... radiative conditions on the total part + relaxation toward climatology
665         ! ... jpjsob,(jpisdp1,jpisfm1)
666         DO jj = njs0, njs1
667            DO jk = 1, jpkm1
668               DO ji = 1, jpi
669                  z05cx = v_cysbnd(ji,jk)
670                  z05cx = z05cx / e2t(ji,jj+1)
671                  z05cx = max( z05cx, -1. )
672         ! ... z05c > 0, inflow  zin=0, ztau=1
673         !         =< 0, outflow zin=1, ztau=rtaus
674                  zin = sign( 1., -1. * z05cx )
675                  zin = 0.5*( zin + abs(zin) )
676                  ztau = (1.-zin )*rtausin + zin * rtaus
677                  z05cx = z05cx * zin
678         ! ... update va with radiative or climatological velocity
679                  va(ji,jj,jk) = va(ji,jj,jk) * (1.-vsmsk(ji,jk)) +             &
680                                 vsmsk(ji,jk) * ( ( 1. + z05cx - ztau )         &
681                                 * vsbnd(ji,jk,nib ,nitm) - 2.*z05cx              &
682                                 * vsbnd(ji,jk,nibm,nit ) + ztau * vfos (ji,jk) ) &
683                                 / (1. - z05cx)
684               END DO
685            END DO
686         END DO
687      END IF
688
689   END SUBROUTINE obc_dyn_south
690#else
691   !!=================================================================================
692   !!                       ***  MODULE  obcdyn  ***
693   !! Ocean dynamics:   Radiation of velocities on each open boundary
694   !!=================================================================================
695CONTAINS
696
697   SUBROUTINE obc_dyn
698                              ! No open boundaries ==> empty routine
699   END SUBROUTINE obc_dyn
700#endif
701
702END MODULE obcdyn
Note: See TracBrowser for help on using the repository browser.