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/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/OBC – NEMO

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications

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.