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/2012/dev_CMCC_INGV_2012/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/2012/dev_CMCC_INGV_2012/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90 @ 3646

Last change on this file since 3646 was 3646, checked in by vichi, 11 years ago

Add the resulting merged branch from CMCC and INGV 2012 developments

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