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/TAM_V3_0/NEMO/OPA_SRC/OBC – NEMO

source: branches/TAM_V3_0/NEMO/OPA_SRC/OBC/obcdyn.F90 @ 1884

Last change on this file since 1884 was 1884, checked in by rblod, 14 years ago

Light adaptation of NEMO direct model routine to handle TAM

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