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.
obcspg.F90 in trunk/NEMO/OPA_SRC/OBC – NEMO

source: trunk/NEMO/OPA_SRC/OBC/obcspg.F90 @ 3

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

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.4 KB
Line 
1MODULE obcspg 
2   !!======================================================================
3   !!                       ***  MODULE  obcspg  ***
4   !! Open Boundaries  :   Radiation of barotropic stream function on each
5   !!                      open boundary
6   !!======================================================================
7#if defined key_obc && defined key_dynspg_rl
8   !!----------------------------------------------------------------------
9   !!   'key_obc'    and                            Open Boundary Condition
10   !!   'key_dynspg_rl'                                 Rigid-Lid
11   !!----------------------------------------------------------------------
12   !!   obc_spg       : call the subroutine for each open boundary
13   !!   obc_spg_east  : radiation of the east open boundary streamfunction
14   !!   obc_spg_west  : radiation of the west open boundary streamfunction
15   !!   obc_spg_north : radiation of the north open boundary streamfunction
16   !!   obc_spg_south : radiation of the south open boundary streamfunction
17   !!----------------------------------------------------------------------
18   !! * Modules used
19   USE oce             ! ocean dynamics and tracers variables
20   USE dom_oce         ! ocean space and time domain variables
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE phycst          ! physical constants
23   USE obc_oce         ! ocean open boundary conditions
24   USE lib_mpp         ! for mppobc
25   USE in_out_manager  ! I/O manager
26
27   IMPLICIT NONE
28   PRIVATE
29
30   !! * Accessibility
31   PUBLIC obc_spg     ! routine called in step.F90 (rigid lid case)
32
33   !! * Module variables
34   INTEGER ::   ji, jj, jk, jnic   ! 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   !! * Substitutions
49#  include "obc_vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !!   OPA 9.0 , LODYC-IPSL  (2003)
52   !!----------------------------------------------------------------------
53
54CONTAINS
55
56   SUBROUTINE obc_spg ( kt )
57      !!----------------------------------------------------------------------
58      !!                    ***  ROUTINE obc_spg  ***
59      !!
60      !! **  Purpose :
61      !!       Compute now barotropic stream function at the open boundaries.
62      !!       (lpeastobc, and/or lpwestobc, and/or lpnorthobc, and/or lpsouthobc).
63      !!       Deduce the correct bsf trend on the open boundaries and isolated
64      !!       coastlines previous to the call to the barotropic solver.
65      !!
66      !! ** Method :
67      !!      In case of open boundaries, there can be a net barotropic flow
68      !!      through the boundaries, hence the potential on the coastlines
69      !!      on each side of the OBC is different.
70      !!      This routine:
71      !!           1. compute the contribution of the isolated coastlines to the
72      !!              rhs of the barotropic equation
73      !!           2. compute the contribution of the OBC to the rhs of the
74      !!              barotropic equation using a radiation equation as explained
75      !!              in the OBC routine.
76      !!
77      !! Reference :
78      !!   Marchesiello P., 1995, these de l'universite J. Fourier, Grenoble, France.
79      !! History :
80      !!        ! 95-03 (J.-M. Molines) Original, SPEM
81      !!        ! 97-07 (G. Madec, J.-M. Molines) additions
82      !!   8.5  ! 02-10 (C. Talandier, A-M. Treguier) F90
83      !!----------------------------------------------------------------------
84      !! * Arguments
85      INTEGER, INTENT( in ) ::   kt
86      !!----------------------------------------------------------------------
87
88      ! 0. Local constant initialization
89      ! --------------------------------
90
91      IF( kt == nit000 .OR. ln_rstart ) THEN
92         ! ... Boundary restoring coefficient
93         rtaue = 2. * rdt / rdpeob
94         rtauw = 2. * rdt / rdpwob
95         rtaun = 2. * rdt / rdpnob
96         rtaus = 2. * rdt / rdpsob
97         ! ... Boundary restoring coefficient for inflow ( all boundaries)
98         rtauein = 2. * rdt / rdpein 
99         rtauwin = 2. * rdt / rdpwin
100         rtaunin = 2. * rdt / rdpnin
101         rtausin = 2. * rdt / rdpsin 
102      END IF
103
104      ! ... right hand side of the barotropic elliptic equation
105      gcbob(:,:) = 0.e0
106
107      ! 1. Isolated coastline contribution to the RHS of the barotropic Eq.
108      ! -------------------------------------------------------------------
109      DO jnic = 1, nbobc-1
110         DO jj = 1, jpj
111            DO ji = 1, jpi
112               gcbob(ji,jj) = gcbob(ji,jj) + gcfobc(ji,jj,jnic) * gcbic(jnic)
113            END DO
114         END DO
115      END DO
116
117      ! 2. East open boundary
118      ! ---------------------
119
120      IF( lpeastobc ) THEN
121         CALL obc_spg_east( kt )
122      END IF
123
124      ! 3. West open boundary
125      ! ---------------------
126
127      IF( lpwestobc ) THEN
128         CALL obc_spg_west( kt )
129      END IF
130       
131      ! 4. North open boundary
132      ! ----------------------
133
134      IF( lpnorthobc ) THEN
135         CALL obc_spg_north( kt )
136      END IF
137 
138      ! 5. South open boundary
139      ! ----------------------
140
141      IF( lpsouthobc ) THEN
142         CALL obc_spg_south( kt )
143      END IF
144 
145# if defined key_mpp
146      CALL mpp_lnk_2d( gcbob, 'G', 1. )
147# endif
148 
149   END SUBROUTINE obc_spg
150
151   SUBROUTINE obc_spg_east ( kt )
152      !!------------------------------------------------------------------------------
153      !!                     SUBROUTINE obc_spg_east
154      !!                    *************************
155      !! ** Purpose :
156      !!      Apply the radiation algorithm on east OBC stream function.
157      !!      If the logical lfbceast is .TRUE., there is no radiation but only fixed OBC
158      !!
159      !!  History :
160      !!         ! 95-03 (J.-M. Molines) Original from SPEM
161      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
162      !!         ! 97-12 (M. Imbard) Mpp adaptation
163      !!         ! 00-06 (J.-M. Molines)
164      !!    8.5  ! 02-10 (C. Talandier, A-M Treguier) F90
165      !!------------------------------------------------------------------------------
166      !! * Arguments
167      INTEGER, INTENT( in ) ::   kt
168
169      !! * Local declarations
170      INTEGER ::   ij
171
172      REAL(wp) ::   z2dtr, ztau, zin
173      REAL(wp) ::   z05cx, zdt, z4nor2, z2dx, z2dy
174
175      !!------------------------------------------------------------------------------
176      !!  OPA 8.5, LODYC-IPSL (2002)
177      !!------------------------------------------------------------------------------
178
179      ! 1. First three time steps and more if lfbceast is .TRUE.
180      !    In that case open boundary conditions are FIXED.
181      ! --------------------------------------------------------
182
183      IF( ( kt < nit000 + 3 .AND. .NOT.ln_rstart) .OR. lfbceast ) THEN
184
185         ! 1.1 Fixed barotropic stream function
186         ! ------------------------------------
187         DO jj = nje0m1, nje1 
188            ij = jj -1 + njmpp
189            bsfeob(jj)=bfoe(ij)
190         END DO
191
192      ELSE
193
194      ! 2. Beyond the fourth time step if lfbceast is .FALSE.
195      ! -----------------------------------------------------
196
197         ! 2.1. Barotropic stream function radiation
198         ! ----------------------------------------
199         !
200         !          nibm2      nibm      nib
201         !            |   nibm  |   nib   |///
202         !            |    |    |    |    |///
203         !  jj-line --f----v----f----v----f---
204         !            |         |         |///
205         !            |    |    |    |    |///
206         !         jpieob-2   jpieob-1    jpieob
207         !                 |         |       
208         !              jpieob-1    jpieob     
209         !
210         ! ... radiative conditions plus restoring term toward climatology
211         ! ... Initialize bsfeob to clim in any case, at the first step
212         !     to ensure proper values at the ends of the open line.
213         ! ... Take care that the j indices starts at nje0 (jpjed) and finish
214         !     at nje1 (jpjef) to be sure that jpjefm1 and jpjef are set OK.
215         DO ji = fs_nie0-1, fs_nie1-1 ! Vector opt.
216            DO jj = nje0p1, nje1m2 
217               ij = jj -1 + njmpp
218         ! ... 2* gradi(bsf) (v-point i=nibm, time mean)
219               z2dx = ( bebnd(jj,nibm ,nit) + bebnd(jj,nibm ,nitm2) - 2.*bebnd(jj,nibm2,nitm) ) &
220                      / e1v(ji,jj)
221         ! ... 2* gradj(bsf) (f-point i=nibm, time nitm)
222               z2dy = ( bebnd(jj+1,nibm,nitm) - bebnd(jj-1,nibm,nitm) ) / e2f(ji,jj)
223         ! ... square of the norm of grad(bsf)
224               z4nor2 = z2dx * z2dx + z2dy * z2dy
225         ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
226               zdt = bebnd(jj,nibm,nitm2) - bebnd(jj,nibm,nit)
227         ! ... i-phase speed ratio (bounded by 1) and MASKED!
228               IF( z4nor2 == 0 ) THEN
229                  IF(lwp) WRITE(numout,*)' PB dans obc_spg_east au pt ',jj,' : z4nor=0'
230                  z4nor2 = 0.001
231               END IF
232               z05cx = zdt * z2dx / z4nor2 * bmask(ji,jj)
233               z05cx = z05cx / e1v(ji+1,jj)
234               z05cx = min( z05cx, 1. )
235         ! ... z05cx < 0, inflow  zin=0, ztau=1 
236         !          => 0, outflow zin=1, ztau=rtaue
237               zin = sign( 1., z05cx )
238               zin = 0.5*( zin + abs(zin) )
239         ! ... Modification JM:  We maintain a restoring term toward
240         !                   bsfb even in case of inflow
241         ! But restoring is stronger when in flow (10 days) (ztau in set in obcspg.F)
242               ztau = (1.-zin ) * rtauein + zin * rtaue
243               z05cx = z05cx * zin
244         ! ... update bsfn with radiative or climatological bsf (not mask!)
245               bsfeob(jj) = ( ( 1. - z05cx - ztau ) * bebnd(jj,nib ,nitm) + 2.*z05cx  &
246                               * bebnd(jj,nibm,nit) + ztau * bfoe (ij) )              &
247                            / (1. + z05cx)
248            END DO
249         END DO
250
251      END IF
252# if defined key_mpp
253      CALL mppobc(bsfeob,jpjed,jpjef,jpieob-1,1,2,jpj)
254# endif
255
256      ! 3. right hand side of the barotropic elliptic equation
257      ! ------------------------------------------------------
258 
259      IF( ( neuler == 0 ) .AND. ( kt == nit000 ) ) THEN
260         z2dtr=1./rdt
261      ELSE
262         z2dtr=1./2./rdt
263      END IF
264      DO ji = fs_nie0-1, fs_nie1-1 ! Vector opt.
265         DO jj = nje0m1, nje1 
266            gcbob(ji,jj) = gcbob(ji,jj) - hvr(ji+1,jj) * e2v(ji+1,jj) / e1v(ji+1,jj) &
267                           * ( bsfeob(jj) - bsfb(ji+1,jj) ) * z2dtr * bmask(ji,jj) 
268         END DO
269      END DO
270
271   END SUBROUTINE obc_spg_east
272
273   SUBROUTINE obc_spg_west ( kt )
274      !!------------------------------------------------------------------------------
275      !!                     SUBROUTINE obc_spg_west
276      !!                    *************************
277      !! ** Purpose :
278      !!      Apply the radiation algorithm on west OBC stream function.
279      !!      If the logical lfbcwest is .TRUE., there is no radiation but only fixed OBC
280      !!
281      !!  History :
282      !!         ! 95-03 (J.-M. Molines) Original from SPEM
283      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
284      !!         ! 97-12 (M. Imbard) Mpp adaptation
285      !!         ! 00-06 (J.-M. Molines)
286      !!    8.5  ! 02-10 (C. Talandier, A-M Treguier) F90
287      !!------------------------------------------------------------------------------
288      !! * Arguments
289      INTEGER, INTENT( in ) ::   kt
290
291      !! * Local declarations
292      INTEGER ::   ij
293
294      REAL(wp) ::   z2dtr, ztau, zin
295      REAL(wp) ::   z05cx, zdt, z4nor2, z2dx, z2dy
296
297      !!------------------------------------------------------------------------------
298      !!  OPA 8.5, LODYC-IPSL (2002)
299      !!------------------------------------------------------------------------------
300
301      ! 1. First three time steps and more if lfbcwest is .TRUE.
302      !    In that case open boundary conditions are FIXED.
303      ! --------------------------------------------------------
304
305      IF( ( kt < nit000 + 3 .AND. .NOT.ln_rstart ) .OR. lfbcwest ) THEN
306
307         ! 1.1 Fixed barotropic stream function
308         ! ------------------------------------
309         DO jj = njw0m1, njw1
310            ij = jj -1 + njmpp
311            bsfwob(jj)=bfow(ij)
312         END DO
313
314      ELSE
315
316      ! 2. Beyond the fourth time step if lfbcwest is .FALSE.
317      ! -----------------------------------------------------
318
319         ! 2.1. Barotropic stream function radiation
320         ! ----------------------------------------
321         !
322         !         nib       nibm     nibm2
323         !       ///|   nib   |   nibm  |
324         !       ///|    |    |    |    |
325         !       ---f----v----f----v----f-- jj-line
326         !       ///|         |         |
327         !       ///|    |    |    |    |
328         !        jpiwob    jpiwob+1    jpiwob+2
329         !               |         |       
330         !             jpiwob+1    jpiwob+2     
331         !
332         ! ... radiative conditions plus restoring term toward climatology
333         ! ... Initialize bsfwob to clim in any case, at the first step
334         !     to ensure proper values at the ends of the open line.
335         ! ... Take care that the j indices starts at njw0 (jpjwd) and finish
336         ! ... at njw1 (jpjwf) to be sure that jpjwfm1 and jpjwf are set OK.
337         DO ji = fs_niw0+1, fs_niw1+1 ! Vector opt.
338            DO jj = njw0p1, njw1m2
339               ij = jj -1 + njmpp
340         ! ... 2* gradi(bsf) (v-point i=nibm, time mean)
341               z2dx = ( -  bwbnd(jj,nibm ,nit ) - bwbnd(jj,nibm ,nitm2) + 2.*bwbnd(jj,nibm2,nitm ) ) &
342                      / e1v(ji+1,jj)
343         ! ... 2* gradj(bsf) (f-point i=nibm, time nitm)
344               z2dy = ( bwbnd(jj+1,nibm,nitm) - bwbnd(jj-1,nibm,nitm) ) / e2f(ji,jj)
345         ! ... square of the norm of grad(bsf)
346               z4nor2 = z2dx * z2dx + z2dy * z2dy
347         ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
348               zdt = bwbnd(jj,nibm,nitm2) - bwbnd(jj,nibm,nit)
349         ! ... i-phase speed ratio (bounded by 1) and MASKED!
350               IF( z4nor2 == 0 ) THEN
351                  IF(lwp) WRITE(numout,*)' PB dans obc_spg_west au pt ',jj,' : z4nor =0'
352                  z4nor2=0.0001
353               END IF
354               z05cx = zdt * z2dx / z4nor2 * bmask(ji,jj)
355               z05cx = z05cx / e1v(ji,jj)
356               z05cx = max( z05cx, -1. )
357         ! ... z05cx => 0, inflow  zin=0, ztau=1 
358         !           <  0, outflow zin=1, ztau=rtauw
359               zin = sign( 1., -1. * z05cx )
360               zin = 0.5*( zin + abs(zin) )
361               ztau = (1.-zin )*rtauwin + zin * rtauw
362               z05cx = z05cx * zin
363         !  ... update bsfn with radiative or climatological bsf (not mask!)
364               bsfwob(jj) = ( ( 1. + z05cx - ztau ) * bwbnd(jj,nib ,nitm) - 2.*z05cx &
365                               * bwbnd(jj,nibm,nit) + ztau * bfow (ij) )             &
366                            / (1. - z05cx)
367            END DO
368         END DO
369
370      END IF
371# if defined key_mpp 
372      CALL mppobc(bsfwob,jpjwd,jpjwf,jpiwob+1,1,2,jpj) 
373# endif 
374
375      ! 3. right hand side of the barotropic elliptic equation
376      ! -------------------------------------------------------
377
378      IF( ( neuler == 0 ) .AND. ( kt == nit000 ) ) THEN
379         z2dtr=1./rdt
380      ELSE
381         z2dtr=1./2./rdt
382      END IF
383      DO ji = fs_niw0+1, fs_niw1+1 ! Vector opt.
384         DO jj = njw0m1, njw1
385            gcbob(ji,jj) = gcbob(ji,jj) - hvr(ji,jj) * e2v(ji,jj) / e1v(ji,jj)     &
386                           * ( bsfwob(jj) - bsfb(ji-1,jj) ) * z2dtr * bmask(ji,jj)
387         END DO
388      END DO
389
390   END SUBROUTINE obc_spg_west
391
392   SUBROUTINE obc_spg_north ( kt )
393      !!------------------------------------------------------------------------------
394      !!                     SUBROUTINE obc_spg_north
395      !!                    *************************
396      !! ** Purpose :
397      !!      Apply the radiation algorithm on north OBC stream function.
398      !!      If the logical lfbcnorth is .TRUE., there is no radiation but only fixed OBC
399      !!
400      !!  History :
401      !!         ! 95-03 (J.-M. Molines) Original from SPEM
402      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
403      !!         ! 97-12 (M. Imbard) Mpp adaptation
404      !!         ! 00-06 (J.-M. Molines)
405      !!    8.5  ! 02-10 (C. Talandier, A-M Treguier) F90
406      !!------------------------------------------------------------------------------
407      !! * Arguments
408      INTEGER, INTENT( in ) ::   kt
409
410      !! * Local declarations
411      INTEGER ::   ii
412
413      REAL(wp) ::   z2dtr, ztau, zin
414      REAL(wp) ::   z05cx, zdt, z4nor2, z2dx, z2dy
415
416      !!------------------------------------------------------------------------------
417      !!  OPA 8.5, LODYC-IPSL (2002)
418      !!------------------------------------------------------------------------------
419
420      ! 1. First three time steps and more if lfbcnorth is .TRUE.
421      !    In that case open boundary conditions are FIXED.
422      ! --------------------------------------------------------
423
424      IF( ( kt < nit000 + 3 .AND. .NOT.ln_rstart ) .OR. lfbcnorth ) THEN
425
426         ! 1.1 Fixed barotropic stream function
427         ! ------------------------------------
428         DO ji = nin0m1, nin1
429            ii = ji -1 + nimpp
430            bsfnob(ji)=bfon(ii)
431         END DO
432
433      ELSE     
434
435      ! 2. Beyond the fourth time step if lfbcnorth is .FALSE.
436      ! -----------------------------------------------------
437
438         ! 2.1. Barotropic stream function radiation
439         ! -----------------------------------------
440         !
441         !           ji-row
442         !             |
443         !        ////////////
444         !        ////////////
445         !   nib  -----f------  jpjnob
446         !             |   
447         !      nib--  u   ---- jpjnob
448         !             |       
449         !  nibm  -----f-----   jpjnob-1
450         !             |       
451         !     nibm--  u   ---- jpjnob-1
452         !             |       
453         !  nibm2 -----f-----   jpjnob-2
454         !             |
455         !
456         ! ... radiative conditions plus restoring term toward climatology
457         ! ... z05cx is always the cross boundary phase velocity
458         ! ... Initialize bsfnob to clim in any case, at the first step
459         !     to ensure proper values at the ends of the open line.
460         ! ... Take care that the i indices starts at nin0 (jpind) and finish
461         ! ... at nin1 (jpinf) to be sure that jpinfm1 and jpinf are set OK.
462         DO jj = fs_njn0-1, fs_njn1-1 ! Vector opt.
463            DO ji = nin0p1, nin1m2
464               ii = ji -1 + nimpp
465         ! ... 2* gradj(bsf) (u-point i=nibm, time mean)
466               z2dx = ( bnbnd(ji,nibm ,nit) + bnbnd(ji,nibm ,nitm2) - 2.*bnbnd(ji,nibm2,nitm) ) &
467                      / e2u(ji,jj)
468         ! ... 2* gradi(bsf) (f-point i=nibm, time nitm)
469               z2dy = ( bnbnd(ji+1,nibm,nitm) - bnbnd(ji-1,nibm,nitm) ) / e1f(ji,jj)
470         ! ... square of the norm of grad(bsf)
471               z4nor2 = z2dx * z2dx + z2dy * z2dy
472         ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
473               zdt = bnbnd(ji,nibm,nitm2) - bnbnd(ji,nibm,nit)
474         ! ... j-phase speed ratio (bounded by 1) and MASKED!
475               IF( z4nor2 == 0 ) THEN
476                  IF(lwp) WRITE(numout,*)' PB dans obc_spg_north au pt',ji,' : z4nor =0'
477               END IF
478               z05cx = zdt * z2dx / z4nor2 * bmask(ji,jj)
479               z05cx = z05cx / e2u(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               ztau = (1.-zin ) * rtaunin + zin * rtaun
486               z05cx = z05cx * zin
487         ! ... update bsfn with radiative or climatological bsf (not mask!)
488               bsfnob(ji) = ( ( 1. - z05cx - ztau ) * bnbnd(ji,nib ,nitm) + 2.*z05cx  &
489                               * bnbnd(ji,nibm,nit) + ztau * bfon (ii) )              &
490                            / (1. + z05cx)
491            END DO
492         END DO
493
494      END IF
495# if defined key_mpp
496      call mppobc(bsfnob,jpind,jpinf,jpjnob-1,1,1,jpi)
497# endif
498
499      ! 3. right hand side of the barotropic elliptic equation
500      !-------------------------------------------------------
501
502      IF( ( neuler == 0 ) .AND. ( kt == nit000 ) ) THEN
503         z2dtr=1./rdt
504      ELSE
505         z2dtr=1./2./rdt
506      END IF
507      DO jj = fs_njn0-1, fs_njn1-1 ! Vector opt.
508         DO ji = nin0m1, nin1
509            gcbob(ji,jj) = gcbob(ji,jj) - hur(ji,jj+1) *  e1u(ji,jj+1) / e2u(ji,jj+1) &
510                           * ( bsfnob(ji) - bsfb(ji,jj+1) ) * z2dtr * bmask(ji,jj)
511         END DO
512      END DO
513
514   END SUBROUTINE obc_spg_north
515
516   SUBROUTINE obc_spg_south ( kt )
517      !!------------------------------------------------------------------------------
518      !!                     SUBROUTINE obc_spg_south
519      !!                    *************************
520      !! ** Purpose :
521      !!      Apply the radiation algorithm on south OBC stream function.
522      !!      If the logical lfbcsouth is .TRUE., there is no radiation but only fixed OBC
523      !!
524      !!  History :
525      !!         ! 95-03 (J.-M. Molines) Original from SPEM
526      !!         ! 97-07 (G. Madec, J.-M. Molines) additions
527      !!         ! 97-12 (M. Imbard) Mpp adaptation
528      !!         ! 00-06 (J.-M. Molines)
529      !!    8.5  ! 02-10 (C. Talandier, A-M Treguier) F90
530      !!------------------------------------------------------------------------------
531      !! * Arguments
532      INTEGER, INTENT(in) ::   kt  ! ocean time-step index
533
534      !! * Local declarations
535      INTEGER  ::   ii             ! temporary integers
536      REAL(wp) ::   &
537         z2dtr, ztau, zin   ,   &  ! temporary scalars
538         z05cx, zdt , z4nor2,   &  !    "         "
539         z2dx , z2dy               !    "         "
540      !!------------------------------------------------------------------------------
541
542      ! 1. First three time steps and more if lfbcsouth is .TRUE.
543      !    In that case open boundary conditions are FIXED.
544      ! --------------------------------------------------------
545     
546      IF( ( kt < nit000 + 3 .AND. .NOT.ln_rstart ) .OR. lfbcsouth ) THEN
547
548         ! 1.1 Fixed barotropic stream function
549         ! ------------------------------------
550         DO ji = nis0m1, nis1 
551            ii = ji -1 + nimpp 
552            bsfsob(ji)=bfos(ii)
553         END DO
554
555      ELSE
556
557      ! 2. Beyond the fourth time step if lfbcsouth is .FALSE.
558      ! -----------------------------------------------------
559
560         ! 2.1. Barotropic stream function radiation
561         ! -----------------------------------------
562         !
563         !           ji-row
564         !             |
565         ! nibm2  -----f------  jpjsob + 2
566         !             |   
567         !   nibm  --  u  ----- jpjsob + 2
568         !             |       
569         !  nibm  -----f-----   jpjsob + 1
570         !             |       
571         !    nib  --  u  ----- jpjsob + 1
572         !             |       
573         !    nib -----f-----   jpjsob
574         !        ///////////     
575         !        ///////////
576         !
577         ! ... radiative conditions plus restoring term toward climatology
578         ! ... z05cx is always the cross boundary phase velocity
579         ! ... Initialize bsfsob to clim in any case, at the first step
580         !     to ensure proper values at the ends of the open line.
581         ! ... Take care that the i indices starts at nis0 (jpisd) and finish
582         ! ... at nis1 (jpisf) to be sure that jpisfm1 and jpisf are set OK.
583         DO jj = fs_njs0+1, fs_njs1+1 ! Vector opt.
584            DO ji = nis0p1, nis1m2
585               ii = ji -1 + nimpp
586         ! ... 2* gradj(bsf) (u-point i=nibm, time mean)
587               z2dx = ( - bsbnd(ji,nibm ,nit) - bsbnd(ji,nibm ,nitm2) + 2.*bsbnd(ji,nibm2,nitm) ) &
588                      / e2u(ji,jj+1)
589         ! ... 2* gradi(bsf) (f-point i=nibm, time nitm)
590               z2dy = ( bsbnd(ji+1,nibm,nitm) - bsbnd(ji-1,nibm,nitm) ) / e1f(ji,jj)
591         ! ... square of the norm of grad(bsf)
592               z4nor2 = z2dx * z2dx + z2dy * z2dy
593         ! ... minus time derivative (leap-frog) at nibm, without / 2 dt
594               zdt = bsbnd(ji,nibm,nitm2) - bsbnd(ji,nibm,nit)
595         ! ... j-phase speed ratio (bounded by -1) and MASKED!
596               IF( z4nor2 == 0 ) THEN
597                  IF(lwp) WRITE(numout,*)' PB dans obc_spg_south au pt ',ji,' : z4nor =0'
598               END IF
599               z05cx = zdt * z2dx / z4nor2 * bmask(ji,jj)
600               z05cx = z05cx / e2u(ji,jj)
601               z05cx = max( z05cx, -1. )
602         ! ... z05cx => 0, inflow  zin=0, ztau=1
603         !           <  0, outflow zin=1, ztau=rtaus
604               zin = sign( 1., -1. * z05cx )
605               zin = 0.5*( zin + abs(zin) )
606               ztau = (1.-zin ) *rtausin  + zin * rtaus
607               z05cx = z05cx * zin
608         ! ... update bsfn with radiative or climatological bsf (not mask!)
609               bsfsob(ji) = ( ( 1. + z05cx - ztau ) * bsbnd(ji,nib ,nitm) - 2.*z05cx  & 
610                               * bsbnd(ji,nibm,nit) + ztau * bfos (ii) )              &
611                            / (1. - z05cx)
612            END DO
613         END DO
614
615      END IF
616# if defined key_mpp
617      CALL mppobc(bsfsob,jpisd,jpisf,jpjsob+1,1,1,jpi)
618# endif
619 
620      ! 3. right hand side of the barotropic elliptic equation
621      ! -------------------------------------------------------
622
623      IF( ( neuler == 0 ) .and. ( kt == nit000 ) ) THEN
624         z2dtr=1./rdt
625      ELSE
626         z2dtr=1./2./rdt
627      END IF
628      DO jj = fs_njs0+1, fs_njs1+1 ! Vector opt.
629         DO ji = nis0m1, nis1 
630            gcbob(ji,jj) = gcbob(ji,jj) - hur(ji,jj) * e1u(ji,jj) / e2u(ji,jj) &
631                           * ( bsfsob(ji) - bsfb(ji,jj-1) ) * z2dtr * bmask(ji,jj)
632         END DO
633      END DO
634
635   END SUBROUTINE obc_spg_south
636
637#else
638   !!----------------------------------------------------------------------
639   !!   Default case :                                         Empty module
640   !!----------------------------------------------------------------------
641CONTAINS
642   SUBROUTINE obc_spg( kt )        ! Empty routine
643      INTEGER, INTENT( in ) :: kt
644      WRITE(*,*) kt
645   END SUBROUTINE obc_spg
646#endif
647
648   !!======================================================================
649END MODULE obcspg
Note: See TracBrowser for help on using the repository browser.