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

source: trunk/NEMO/OPA_SRC/OBC/obcini.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 31.1 KB
Line 
1 MODULE obcini
2   !!=================================================================================
3   !!                       ***  MODULE  obcini  ***
4   !! OBC initial state :  Open boundary initial state
5   !!=================================================================================
6#if defined key_obc
7   !!---------------------------------------------------------------------------------
8   !!   'key_obc'                                             Open Boundary Conditions
9   !!---------------------------------------------------------------------------------
10   !!   obc_init       : initialization for the open boundary condition
11   !!---------------------------------------------------------------------------------
12   !! * Modules used
13   USE oce             ! ocean dynamics and tracers variables
14   USE dom_oce         ! ocean space and time domain variables
15   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
16   USE phycst          ! physical constants
17   USE obc_oce         ! ocean open boundary conditions
18   USE lib_mpp         ! for mpp_sum
19   USE in_out_manager  ! I/O units
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! * Routine accessibility
25   PUBLIC obc_init        ! routine called by opa.F90
26
27   !! * Substitutions
28#  include "obc_vectopt_loop_substitute.h90"
29   !!---------------------------------------------------------------------------------
30   !!   OPA 9.0 , LOCEAN-IPSL (2005)
31   !! $Header$
32   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
33   !!---------------------------------------------------------------------------------
34
35CONTAINS
36   
37   SUBROUTINE obc_init
38      !!----------------------------------------------------------------------
39      !!                 ***  ROUTINE obc_init  ***
40      !!         
41      !! ** Purpose :   Initialization of the dynamics and tracer fields at
42      !!      the open boundaries.
43      !!
44      !! ** Method  :   initialization of open boundary variables
45      !!      (u, v, bsf) over 3 time step and 3 rows
46      !!      (t, s)      over 2 time step and 2 rows
47      !!      if ln_rstart = .FALSE. : no restart, fields set to zero
48      !!      if ln_rstart = .TRUE.  : restart, fields are read in a file
49      !!      if rdpxxx = 0 then lfbc is set true for this boundary.
50      !!
51      !! ** Input   :   restart.obc file, restart file for open boundaries
52      !!
53      !! History :
54      !!   8.0  !  97-07  (G. Madec)  Original code
55      !!        !  97-11  (J.M. Molines)
56      !!   8.5  !  02-11  (C. Talandier, A-M. Treguier) Free surface, F90
57      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
58      !!----------------------------------------------------------------------
59      !! * Modules used
60      USE obcrst,   ONLY :   obc_rst_lec   ! Make obc_rst_lec routine available
61      USE obcdom,   ONLY :   obc_dom       ! Make obc_dom routine available
62
63      !! * Local declarations
64      INTEGER  ::   ji, jj, istop , inumfbc
65      INTEGER, DIMENSION(4) ::   icorner
66      REAL(wp) ::   zbsic1, zbsic2, zbsic3
67      REAL(wp), DIMENSION(2) ::   ztestmask
68
69      NAMELIST/namobc/ rdpein, rdpwin, rdpnin, rdpsin,   &
70         &             rdpeob, rdpwob, rdpnob, rdpsob,   &
71         &             zbsic1, zbsic2, zbsic3,           &
72         &             nbic, volemp, nobc_dta,           &
73         &             ln_obc_clim, ln_vol_cst, ln_obc_fla
74      !!----------------------------------------------------------------------
75
76      IF(lwp) WRITE(numout,*)
77      IF(lwp) WRITE(numout,*) 'obc_init : initialization of open boundaries'
78      IF(lwp) WRITE(numout,*) '~~~~~~~~'
79
80
81      ! 0. read namelist parameters
82      ! ---------------------------
83      ! default values already set except:
84      zbsic1 = 0.e0
85      zbsic2 = 0.e0
86      zbsic3 = 0.e0
87
88      ! Namelist namobc : open boundaries
89      REWIND( numnam )
90      READ  ( numnam, namobc )
91
92      bsfic0(1) = zbsic1
93      bsfic (2) = zbsic2
94      bsfic (3) = zbsic3
95
96      ! By security we set rdpxin and rdpxob respectively
97      ! to 1. and 15. if the corresponding OBC is not activated
98      IF( .NOT.lp_obc_east ) THEN
99         rdpein = 1.
100         rdpeob = 15.
101      END IF
102      IF( .NOT.lp_obc_west ) THEN
103         rdpwin = 1.
104         rdpwob = 15.
105      END IF
106      IF( .NOT.lp_obc_north ) THEN
107         rdpnin = 1.
108         rdpnob = 15.
109      END IF
110      IF( .NOT.lp_obc_south ) THEN
111         rdpsin = 1.
112         rdpsob = 15.
113      END IF
114
115      ! number of open boudaries and open boundary indicators
116      nbobc = 0
117      IF( lp_obc_east  )   nbobc = nbobc + 1
118      IF( lp_obc_west  )   nbobc = nbobc + 1
119      IF( lp_obc_north )   nbobc = nbobc + 1
120      IF( lp_obc_south )   nbobc = nbobc + 1
121
122      IF(lwp) WRITE(numout,*) '         Number of open boundaries    nbobc = ',nbobc
123      IF(lwp) WRITE(numout,*)
124      IF( nbobc /= 0 .AND. jperio /= 0 ) &
125           &   CALL ctl_stop( ' Cyclic or symmetric, and open boundary condition are not compatible' )
126
127      ! control prints
128      IF(lwp) WRITE(numout,*) '         namobc'
129      IF(lwp) WRITE(numout,*) ' '
130      IF(lwp) WRITE(numout,*) '         data in file (=1) or     nobc_dta = ', nobc_dta
131      IF(lwp) WRITE(numout,*) '         initial state used (=0)             '
132      IF(lwp) WRITE(numout,*) '         climatology (true) or not:', ln_obc_clim
133      IF(lwp) WRITE(numout,*) ' '
134      IF(lwp) WRITE(numout,*) '                                 WARNING                                                  '
135      IF(lwp) WRITE(numout,*) '         Flather"s algorithm is applied with explicit free surface scheme                 '
136      IF(lwp) WRITE(numout,*) '         or with free surface time-splitting scheme                                       '
137      IF(lwp) WRITE(numout,*) '         Nor radiation neither relaxation is allowed with explicit free surface scheme:   '
138      IF(lwp) WRITE(numout,*) '         Radiation and/or relaxation is allowed with free surface time-splitting scheme '
139      IF(lwp) WRITE(numout,*) '         depending of the choice of rdpXin = rdpXob  = 0. for open boundaries             '
140      IF(lwp) WRITE(numout,*) ' '
141      IF(lwp) WRITE(numout,*) '         For the rigid-lid case or the filtered free surface case,                        '
142      IF(lwp) WRITE(numout,*) '         radiation, relaxation or presciption of data can be applied                      '
143      IF( lwp.AND.lp_obc_east ) THEN
144         WRITE(numout,*) '         East open boundary :'
145         WRITE(numout,*) '              i index                    jpieob = ', jpieob
146         WRITE(numout,*) '              damping time scale (days)  rdpeob = ', rdpeob
147         WRITE(numout,*) '              damping time scale (days)  rdpein = ', rdpein
148      ENDIF
149
150      IF( lwp.AND.lp_obc_west ) THEN
151         WRITE(numout,*) '         West open boundary :'
152         WRITE(numout,*) '              i index                    jpiwob = ', jpiwob
153         WRITE(numout,*) '              damping time scale (days)  rdpwob = ', rdpwob
154         WRITE(numout,*) '              damping time scale (days)  rdpwin = ', rdpwin
155      ENDIF
156
157      IF( lwp.AND.lp_obc_north ) THEN
158         WRITE(numout,*) '         North open boundary :'
159         WRITE(numout,*) '               j index                    jpjnob = ', jpjnob
160         WRITE(numout,*) '               damping time scale (days)  rdpnob = ', rdpnob
161         WRITE(numout,*) '               damping time scale (days)  rdpnin = ', rdpnin
162      ENDIF
163
164      IF( lwp.AND.lp_obc_south ) THEN
165         WRITE(numout,*) '         South open boundary :'
166         WRITE(numout,*) '               j index                    jpjsob = ', jpjsob
167         WRITE(numout,*) '               damping time scale (days)  rdpsob = ', rdpsob
168         WRITE(numout,*) '               damping time scale (days)  rdpsin = ', rdpsin
169         WRITE(numout,*) ' '
170      ENDIF
171
172      ! 1. Initialisation of constants
173      ! ------------------------------
174
175      ! ... convert rdp$ob in seconds
176      rdpein = rdpein * rday
177      rdpwin = rdpwin * rday
178      rdpnin = rdpnin * rday
179      rdpsin = rdpsin * rday
180      rdpeob = rdpeob * rday
181      rdpwob = rdpwob * rday
182      rdpnob = rdpnob * rday
183      rdpsob = rdpsob * rday
184      lfbceast  = .FALSE.
185      lfbcwest  = .FALSE.
186      lfbcnorth = .FALSE.
187      lfbcsouth = .FALSE.
188      inumfbc = 0
189      ! ... look for Fixed Boundaries (rdp = 0 )
190      ! ... When specified, lbcxxx flags are set to TRUE and rdpxxx are set to
191      ! ...  a small arbitrary value, (to avoid division by zero further on).
192      ! ...  rdpxxx is not used anymore.
193      IF( lp_obc_east )  THEN
194         IF( (rdpein+rdpeob) == 0 )  THEN
195            lfbceast = .TRUE.
196            rdpein = 1e-3
197            rdpeob = 1e-3
198            inumfbc = inumfbc+1
199         ELSEIF ( (rdpein*rdpeob) == 0 )  THEN
200            CALL ctl_stop( 'obc_init : rdpein & rdpeob must be both zero or non zero' )
201         END IF
202      END IF
203      IF( lp_obc_west )  THEN
204         IF( (rdpwin + rdpwob) == 0 )  THEN
205            lfbcwest = .TRUE.
206            rdpwin = 1e-3
207            rdpwob = 1e-3
208            inumfbc = inumfbc+1
209         ELSEIF ( (rdpwin*rdpwob) == 0 )  THEN
210            CALL ctl_stop( 'obc_init : rdpwin & rdpwob must be both zero or non zero' )
211         END IF
212      END IF
213      IF( lp_obc_north )  THEN
214         IF( (rdpnin + rdpnob) == 0 )  THEN
215            lfbcnorth = .TRUE.
216            rdpnin = 1e-3
217            rdpnob = 1e-3
218            inumfbc = inumfbc+1
219         ELSEIF ( (rdpnin*rdpnob) == 0 )  THEN
220            CALL ctl_stop( 'obc_init : rdpnin & rdpnob must be both zero or non zero' )
221         END IF
222      END IF
223      IF( lp_obc_south )  THEN
224         IF( (rdpsin + rdpsob) == 0 )  THEN
225            lfbcsouth = .TRUE.
226            rdpsin = 1e-3
227            rdpsob = 1e-3
228            inumfbc = inumfbc+1
229         ELSEIF ( (rdpsin*rdpsob) == 0 )  THEN
230            CALL ctl_stop( 'obc_init : rdpsin & rdpsob must be both zero or non zero' )
231         END IF
232      END IF
233
234      ! 2.  Clever mpp indices for loops on the open boundaries.
235      !     The loops will be performed only on the processors
236      !     that contain a given open boundary.
237      ! --------------------------------------------------------
238
239      IF( lp_obc_east ) THEN
240         ! ...   mpp initialization
241         nie0   = max( 1, min(jpieob   - nimpp+1, jpi     ) )
242         nie1   = max( 0, min(jpieob   - nimpp+1, jpi - 1 ) )
243         nie0p1 = max( 1, min(jpieob+1 - nimpp+1, jpi     ) )
244         nie1p1 = max( 0, min(jpieob+1 - nimpp+1, jpi - 1 ) )
245         nie0m1 = max( 1, min(jpieob-1 - nimpp+1, jpi     ) )
246         nie1m1 = max( 0, min(jpieob-1 - nimpp+1, jpi - 1 ) )
247         nje0   = max( 2, min(jpjed    - njmpp+1, jpj     ) )
248         nje1   = max( 0, min(jpjef    - njmpp+1, jpj - 1 ) )
249         nje0p1 = max( 1, min(jpjedp1  - njmpp+1, jpj     ) )
250         nje0m1 = max( 1, min(jpjed    - njmpp+1, jpj     ) )
251         nje1m1 = max( 0, min(jpjefm1  - njmpp+1, jpj - 1 ) )
252         nje1m2 = max( 0, min(jpjefm1-1- njmpp+1, jpj - 1 ) )
253         IF(lwp) THEN
254            IF( lfbceast ) THEN
255               WRITE(numout,*)'     '
256               WRITE(numout,*)'         Specified East Open Boundary'
257            ELSE
258               WRITE(numout,*)'     '
259               WRITE(numout,*)'         Radiative East Open Boundary'
260            END IF
261         END IF
262      END IF
263
264      IF( lp_obc_west ) THEN
265         ! ...   mpp initialization
266         niw0   = max( 1, min(jpiwob   - nimpp+1, jpi     ) )
267         niw1   = max( 0, min(jpiwob   - nimpp+1, jpi - 1 ) )
268         niw0p1 = max( 1, min(jpiwob+1 - nimpp+1, jpi     ) )
269         niw1p1 = max( 0, min(jpiwob+1 - nimpp+1, jpi - 1 ) )
270         njw0   = max( 2, min(jpjwd    - njmpp+1, jpj     ) )
271         njw1   = max( 0, min(jpjwf    - njmpp+1, jpj - 1 ) )
272         njw0p1 = max( 1, min(jpjwdp1  - njmpp+1, jpj     ) )
273         njw0m1 = max( 1, min(jpjwd    - njmpp+1, jpj     ) )
274         njw1m1 = max( 0, min(jpjwfm1  - njmpp+1, jpj - 1 ) )
275         njw1m2 = max( 0, min(jpjwfm1-1- njmpp+1, jpj - 1 ) )
276         IF(lwp) THEN
277            IF( lfbcwest ) THEN
278               WRITE(numout,*)'     '
279               WRITE(numout,*)'         Specified West Open Boundary'
280            ELSE
281               WRITE(numout,*)'     '
282               WRITE(numout,*)'         Radiative West Open Boundary'
283            END IF
284         END IF
285      END IF
286 
287      IF( lp_obc_north ) THEN
288         ! ...   mpp initialization
289         nin0   = max( 2, min(jpind    - nimpp+1, jpi     ) )
290         nin1   = max( 0, min(jpinf    - nimpp+1, jpi - 1 ) )
291         nin0p1 = max( 1, min(jpindp1  - nimpp+1, jpi     ) )
292         nin0m1 = max( 1, min(jpind    - nimpp+1, jpi     ) )
293         nin1m1 = max( 0, min(jpinfm1  - nimpp+1, jpi - 1 ) )
294         nin1m2 = max( 0, min(jpinfm1-1- nimpp+1, jpi - 1 ) )
295         njn0   = max( 1, min(jpjnob   - njmpp+1, jpj     ) )
296         njn1   = max( 0, min(jpjnob   - njmpp+1, jpj - 1 ) )
297         njn0p1 = max( 1, min(jpjnob+1 - njmpp+1, jpj     ) )
298         njn1p1 = max( 0, min(jpjnob+1 - njmpp+1, jpj - 1 ) )
299         njn0m1 = max( 1, min(jpjnob-1 - njmpp+1, jpj     ) )
300         njn1m1 = max( 0, min(jpjnob-1 - njmpp+1, jpj - 1 ) )
301         IF(lwp) THEN
302            IF( lfbcnorth ) THEN
303               WRITE(numout,*)'     '
304               WRITE(numout,*)'         Specified North Open Boundary'
305            ELSE
306               WRITE(numout,*)'     '
307               WRITE(numout,*)'         Radiative North Open Boundary'
308            END IF
309         END IF
310      END IF
311
312      IF( lp_obc_south ) THEN
313         ! ...   mpp initialization
314         nis0   = max( 2, min(jpisd    - nimpp+1, jpi     ) )
315         nis1   = max( 0, min(jpisf    - nimpp+1, jpi - 1 ) )
316         nis0p1 = max( 1, min(jpisdp1  - nimpp+1, jpi     ) )
317         nis0m1 = max( 1, min(jpisd    - nimpp+1, jpi     ) )
318         nis1m1 = max( 0, min(jpisfm1  - nimpp+1, jpi - 1 ) )
319         nis1m2 = max( 0, min(jpisfm1-1- nimpp+1, jpi - 1 ) )
320         njs0   = max( 1, min(jpjsob   - njmpp+1, jpj     ) )
321         njs1   = max( 0, min(jpjsob   - njmpp+1, jpj - 1 ) )
322         njs0p1 = max( 1, min(jpjsob+1 - njmpp+1, jpj     ) )
323         njs1p1 = max( 0, min(jpjsob+1 - njmpp+1, jpj - 1 ) )
324         IF(lwp) THEN
325            IF( lfbcsouth ) THEN
326               WRITE(numout,*)'     '
327               WRITE(numout,*)'         Specified South Open Boundary'
328            ELSE
329               WRITE(numout,*)'     '
330               WRITE(numout,*)'         Radiative South Open Boundary'
331            END IF
332         END IF
333      END IF
334
335      ! 3. mask correction for OBCs
336      ! ---------------------------
337
338      IF( lp_obc_east ) THEN
339         !... (jpjed,jpjefm1),jpieob
340         DO jj = nje0, nje1m1
341# if defined key_dynspg_rl
342            DO ji = nie0, nie1
343# else
344            DO ji = nie0p1, nie1p1
345# endif
346               bmask(ji,jj) = 0.e0
347            END DO
348         END DO
349
350         ! ... initilization to zero
351         uemsk(:,:) = 0.e0
352         vemsk(:,:) = 0.e0
353         temsk(:,:) = 0.e0
354
355         ! ... set 2D mask on East OBC,  Vopt
356         DO ji = fs_nie0, fs_nie1
357            DO jj = nje0, nje1
358               uemsk(jj,:) = umask(ji,  jj,:)
359               vemsk(jj,:) = vmask(ji+1,jj,:)
360               temsk(jj,:) = tmask(ji+1,jj,:)
361            END DO
362         END DO
363
364      END IF
365
366      IF( lp_obc_west ) THEN
367         ! ... (jpjwd,jpjwfm1),jpiwob
368         DO jj = njw0, njw1m1
369            DO ji = niw0, niw1
370               bmask(ji,jj) = 0.e0
371            END DO
372         END DO
373
374         ! ... initilization to zero
375         uwmsk(:,:) = 0.e0
376         vwmsk(:,:) = 0.e0
377         twmsk(:,:) = 0.e0
378
379         ! ... set 2D mask on West OBC,  Vopt
380         DO ji = fs_niw0, fs_niw1
381            DO jj = njw0, njw1
382               uwmsk(jj,:) = umask(ji,jj,:)
383               vwmsk(jj,:) = vmask(ji,jj,:)
384               twmsk(jj,:) = tmask(ji,jj,:)
385            END DO
386         END DO
387
388      END IF
389
390      IF( lp_obc_north ) THEN
391         ! ... jpjnob,(jpind,jpisfm1)
392# if defined key_dynspg_rl
393         DO jj = njn0, njn1
394# else
395         DO jj = njn0p1, njn1p1
396# endif
397            DO ji = nin0, nin1m1
398               bmask(ji,jj) = 0.e0
399            END DO
400         END DO
401
402         ! ... initilization to zero
403         unmsk(:,:) = 0.e0
404         vnmsk(:,:) = 0.e0
405         tnmsk(:,:) = 0.e0
406
407         ! ... set 2D mask on North OBC,  Vopt
408         DO jj = fs_njn0, fs_njn1
409            DO ji = nin0, nin1
410               unmsk(ji,:) = umask(ji,jj+1,:)
411               vnmsk(ji,:) = vmask(ji,jj  ,:)
412               tnmsk(ji,:) = tmask(ji,jj+1,:)
413            END DO
414         END DO
415
416      END IF
417
418      IF( lp_obc_south ) THEN 
419         ! ... jpjsob,(jpisd,jpisfm1)
420         DO jj = njs0, njs1
421            DO ji = nis0, nis1m1
422               bmask(ji,jj) = 0.e0
423            END DO
424         END DO
425
426         ! ... initilization to zero
427         usmsk(:,:) = 0.e0
428         vsmsk(:,:) = 0.e0
429         tsmsk(:,:) = 0.e0
430
431         ! ... set 2D mask on South OBC,  Vopt
432         DO jj = njs0, njs1
433            DO ji = nis0, nis1
434               usmsk(ji,:) = umask(ji,jj,:)
435               vsmsk(ji,:) = vmask(ji,jj,:)
436               tsmsk(ji,:) = tmask(ji,jj,:)
437            END DO
438         END DO
439
440      END IF
441
442# if defined key_dynspg_flt
443
444      ! ... Initialize obcumask and obcvmask for the Force filtering
445      !     boundary condition in dynspg_flt
446      obcumask(:,:) = umask(:,:,1)
447      obcvmask(:,:) = vmask(:,:,1)
448
449      ! ... Initialize obctmsk on overlap region and obcs. This mask
450      !     is used in obcvol.F90 to calculate cumulate flux E-P.
451      !     - no flux E-P on obcs and overlap region (jpereci = jprecj = 1)
452      obctmsk(:,:) = tmask(:,:,1)     
453      obctmsk(1  ,:) = 0.e0           
454      obctmsk(jpi,:) = 0.e0           
455      obctmsk(:  ,1) = 0.e0           
456      obctmsk(:,jpj) = 0.e0           
457
458      IF( lp_obc_east ) THEN
459         ! ... East obc Force filtering mask for the grad D
460         DO ji = nie0, nie1
461            DO jj = nje0p1, nje1m1
462               obcumask(ji  ,jj)=0.e0
463               obcvmask(ji+1,jj)=0.e0
464            END DO
465         END DO
466
467         ! ... set to 0 on East OBC
468         DO jj = nje0p1, nje1m1
469            DO ji = nie0p1, nie1p1
470               obctmsk(ji,jj) = 0.e0
471            END DO
472         END DO
473      END IF
474
475      IF( lp_obc_west ) THEN
476         ! ... West obc Force filtering mask for the grad D
477         DO ji = niw0, niw1
478            DO jj = njw0p1, njw1m1
479               obcumask(ji,jj)=0.e0
480               obcvmask(ji,jj)=0.e0
481            END DO
482         END DO
483
484         ! ... set to 0 on West OBC
485         DO jj = njw0p1, njw1m1
486            DO ji = niw0, niw1
487               obctmsk(ji,jj) = 0.e0
488            END DO
489         END DO
490      END IF
491
492      IF( lp_obc_north ) THEN
493         ! ... North obc Force filtering mask for the grad D
494         DO jj = njn0, njn1
495            DO ji = nin0p1, nin1m1
496               obcvmask(ji,jj  )=0.e0
497               obcumask(ji,jj+1)=0.e0
498            END DO
499         END DO
500
501         ! ... set to 0 on North OBC
502         DO jj = njn0p1, njn1p1
503            DO ji = nin0p1, nin1m1
504               obctmsk(ji,jj) = 0.e0
505            END DO
506         END DO
507      END IF
508
509      IF( lp_obc_south ) THEN
510         ! ... South obc Force filtering mask for the grad D
511         DO jj = njs0, njs1 
512            DO ji = nis0p1, nis1m1
513               obcumask(ji,jj)=0.e0
514               obcvmask(ji,jj)=0.e0
515            END DO
516         END DO
517
518         ! ... set to 0 on South OBC
519         DO jj = njs0, njs1
520            DO ji = nis0p1, nis1m1
521               obctmsk(ji,jj) = 0.e0
522            END DO
523         END DO
524      END IF
525
526# endif
527
528# if ! defined key_dynspg_rl
529
530      IF ( ln_vol_cst ) THEN
531
532         ! 3.1 Total lateral surface for each open boundary
533         ! ------------------------------------------------
534
535         ! ... West open boundary surface
536         IF( lp_obc_west ) THEN
537            DO ji = niw0, niw1
538               DO jj = 1, jpj 
539                  obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uwmsk(jj,1)
540               END DO
541            END DO
542         END IF
543
544         ! ... East open boundary surface
545         IF( lp_obc_east ) THEN
546            DO ji = nie0, nie1
547               DO jj = 1, jpj 
548                  obcsurftot = obcsurftot+hu(ji,jj)*e2u(ji,jj)*uemsk(jj,1)
549               END DO
550            END DO
551         END IF
552
553         ! ... North open boundary vertical surface
554         IF( lp_obc_north ) THEN
555            DO jj = njn0, njn1
556               DO ji = 1, jpi
557                  obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vnmsk(ji,1)
558               END DO
559            END DO
560         END IF
561
562         ! ... South open boundary vertical surface
563         IF( lp_obc_south ) THEN
564            DO jj = njs0, njs1
565               DO ji = 1, jpi
566                  obcsurftot = obcsurftot+hv(ji,jj)*e1v(ji,jj)*vsmsk(ji,1)
567               END DO
568            END DO
569         END IF
570         IF( lk_mpp )   CALL mpp_sum( obcsurftot )   ! sum over the global domain
571      ENDIF
572# endif
573
574      ! 5. Control print on mask
575      !    The extremities of the open boundaries must be in land
576      !    or else correspond to an "ocean corner" between two open boundaries.
577      !    corner 1 is southwest, 2 is south east, 3 is northeast, 4 is northwest.
578      ! --------------------------------------------------------------------------
579
580      icorner(:)=0
581
582      ! ... control of the west boundary
583      IF( lp_obc_west ) THEN
584         IF( jpiwob < 2 .OR.  jpiwob >= jpiglo-2 ) THEN
585            WRITE(ctmp1,*) ' jpiwob exceed ', jpiglo-2, 'or less than 2'
586            CALL ctl_stop( ctmp1 )
587         END IF
588         ztestmask(:)=0.
589         DO ji=niw0,niw1
590            IF( (njw0 + njmpp - 1) == jpjwd ) ztestmask(1)=ztestmask(1)+ tmask(ji,njw0,1)
591            IF( (njw1 + njmpp - 1) == jpjwf ) ztestmask(2)=ztestmask(2)+ tmask(ji,njw1,1)
592         END DO
593         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
594
595         IF( ztestmask(1) /= 0. ) icorner(1)=icorner(1)+1
596         IF( ztestmask(2) /= 0. ) icorner(4)=icorner(4)+1
597      END IF
598
599      ! ... control of the east boundary
600      IF( lp_obc_east ) THEN
601         IF( jpieob < 4 .OR.  jpieob >= jpiglo ) THEN
602            WRITE(ctmp1,*) ' jpieob exceed ', jpiglo, ' or less than 4'
603            CALL ctl_stop( ctmp1 )
604         END IF
605         ztestmask(:)=0.
606         DO ji=nie0p1,nie1p1
607            IF( (nje0 + njmpp - 1) == jpjed ) ztestmask(1)=ztestmask(1)+ tmask(ji,nje0,1)
608            IF( (nje1 + njmpp - 1) == jpjef ) ztestmask(2)=ztestmask(2)+ tmask(ji,nje1,1)
609         END DO
610         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
611
612        IF( ztestmask(1) /= 0. ) icorner(2)=icorner(2)+1
613        IF( ztestmask(2) /= 0. ) icorner(3)=icorner(3)+1
614      END IF
615
616      ! ... control of the north boundary
617      IF( lp_obc_north ) THEN
618         IF( jpjnob < 4 .OR.  jpjnob >= jpjglo ) THEN
619            WRITE(ctmp1,*) 'jpjnob exceed ', jpjglo, ' or less than 4'
620            CALL ctl_stop( ctmp1 )
621         END IF
622         ztestmask(:)=0.
623         DO jj=njn0p1,njn1p1
624            IF( (nin0 + nimpp - 1) == jpind ) ztestmask(1)=ztestmask(1)+ tmask(nin0,jj,1)
625            IF( (nin1 + nimpp - 1) == jpinf ) ztestmask(2)=ztestmask(2)+ tmask(nin1,jj,1)
626         END DO
627         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
628
629         IF( ztestmask(1) /= 0. ) icorner(4)=icorner(4)+1
630         IF( ztestmask(2) /= 0. ) icorner(3)=icorner(3)+1
631      END IF
632
633      ! ... control of the south boundary
634      IF( lp_obc_south ) THEN
635         IF( jpjsob < 2 .OR.  jpjsob >= jpjglo-2 ) THEN
636            WRITE(ctmp1,*) ' jpjsob exceed ', jpjglo-2, ' or less than 2'
637            CALL ctl_stop( ctmp1 )
638         END IF
639         ztestmask(:)=0.
640         DO jj=njs0,njs1
641            IF( (nis0 + nimpp - 1) == jpisd ) ztestmask(1)=ztestmask(1)+ tmask(nis0,jj,1)
642            IF( (nis1 + nimpp - 1) == jpisf ) ztestmask(2)=ztestmask(2)+ tmask(nis1,jj,1)
643         END DO
644         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
645
646         IF( ztestmask(1) /= 0. ) icorner(1)=icorner(1)+1
647         IF( ztestmask(2) /= 0. ) icorner(2)=icorner(2)+1
648      END IF
649
650      IF( icorner(1) == 2 ) THEN
651         IF(lwp) WRITE(numout,*)
652         IF(lwp) WRITE(numout,*) ' South West ocean corner, two open boudaries'
653         IF(lwp) WRITE(numout,*) ' ========== '
654         IF(lwp) WRITE(numout,*)
655         IF( jpisd /= jpiwob.OR.jpjsob /= jpjwd ) &
656              &   CALL ctl_stop( ' Open boundaries do not fit, we stop' )
657
658      ELSE IF( icorner(1) == 1 ) THEN
659         CALL ctl_stop( ' Open boundaries do not fit at SW corner, we stop' )
660      END IF
661
662      IF( icorner(2) == 2 ) THEN
663          IF(lwp) WRITE(numout,*)
664          IF(lwp) WRITE(numout,*) ' South East ocean corner, two open boudaries'
665          IF(lwp) WRITE(numout,*) ' ========== '
666          IF(lwp) WRITE(numout,*)
667          IF( jpisf /= jpieob+1.OR.jpjsob /= jpjed ) &
668               &   CALL ctl_stop( ' Open boundaries do not fit, we stop' )
669      ELSE IF( icorner(2) == 1 ) THEN
670         CALL ctl_stop( ' Open boundaries do not fit at SE corner, we stop' )
671      END IF
672
673      IF( icorner(3) == 2 ) THEN
674         IF(lwp) WRITE(numout,*)
675         IF(lwp) WRITE(numout,*) ' North East ocean corner, two open boudaries'
676         IF(lwp) WRITE(numout,*) ' ========== '
677         IF(lwp) WRITE(numout,*)
678         IF( jpinf /= jpieob+1 .OR. jpjnob+1 /= jpjef ) &
679              &   CALL ctl_stop( ' Open boundaries do not fit, we stop' )
680       ELSE IF( icorner(3) == 1 ) THEN
681          CALL ctl_stop( ' Open boundaries do not fit at NE corner, we stop' )
682       END IF
683
684      IF( icorner(4) == 2 ) THEN
685         IF(lwp) WRITE(numout,*)
686         IF(lwp) WRITE(numout,*) ' North West ocean corner, two open boudaries'
687         IF(lwp) WRITE(numout,*) ' ========== '
688         IF(lwp) WRITE(numout,*)
689         IF( jpind /= jpiwob.OR.jpjnob+1 /= jpjwf ) &
690              &   CALL ctl_stop( ' Open boundaries do not fit, we stop' )
691       ELSE IF( icorner(4) == 1 ) THEN
692          CALL ctl_stop( ' Open boundaries do not fit at NW corner, we stop' )
693       END IF 
694
695      ! 6. Initialization of open boundary variables (u, v, bsf, t, s)
696      ! --------------------------------------------------------------
697      !   only if at least one boundary is  radiative
698
699      ! ... Restart from restart.obc
700      IF ( inumfbc < nbobc .AND.  ln_rstart ) THEN
701         CALL obc_rst_lec
702      ELSE
703
704          ! ... Initialization to zero of radiation arrays.
705          !     Those have dimensions of local subdomains
706
707          bebnd(:,:,:)   = 0.e0   ;   bnbnd(:,:,:)   = 0.e0
708          uebnd(:,:,:,:) = 0.e0   ;   unbnd(:,:,:,:) = 0.e0
709          vebnd(:,:,:,:) = 0.e0   ;   vnbnd(:,:,:,:) = 0.e0
710          tebnd(:,:,:,:) = 0.e0   ;   tnbnd(:,:,:,:) = 0.e0 
711          sebnd(:,:,:,:) = 0.e0   ;   snbnd(:,:,:,:) = 0.e0
712
713          bwbnd(:,:,:)   = 0.e0   ;   bsbnd(:,:,:)   = 0.e0
714          uwbnd(:,:,:,:) = 0.e0   ;   usbnd(:,:,:,:) = 0.e0
715          vwbnd(:,:,:,:) = 0.e0   ;   vsbnd(:,:,:,:) = 0.e0
716          twbnd(:,:,:,:) = 0.e0   ;   tsbnd(:,:,:,:) = 0.e0 
717          swbnd(:,:,:,:) = 0.e0   ;   ssbnd(:,:,:,:) = 0.e0
718
719      END IF
720
721# if defined key_dynspg_rl
722      ! 7. Isolated coastline arrays initialization (rigid lid case only)
723      ! -----------------------------------------------------------------
724      CALL obc_dom
725# endif
726
727      ! 8. Control print
728      ! ... control of the east boundary
729      IF( lp_obc_east ) THEN
730         istop = 0
731         IF( jpieob < 4 .OR.  jpieob >= jpiglo ) THEN
732            IF(lwp) WRITE(numout,cform_err)
733            IF(lwp) WRITE(numout,*) '            jpieob exceed ', jpim1, ' or less than 4'
734            istop = istop + 1
735         END IF
736
737         IF( lk_mpp ) THEN
738            ! ...
739            IF( nimpp > jpieob-5) THEN
740               IF(lwp) WRITE(numout,cform_err)
741               IF(lwp) WRITE(numout,*) '        A sub-domain is too close to the East OBC'
742               IF(lwp) WRITE(numout,*) '        nimpp must be < jpieob-5'
743               istop = istop + 1
744            ENDIF
745         ELSE
746
747            ! ... stop if  e r r o r (s)   detected
748            IF( istop /= 0 ) THEN
749               WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop'
750               CALL ctl_stop( ctmp1 )
751            ENDIF
752         ENDIF
753      ENDIF
754
755      ! ... control of the west boundary
756      IF( lp_obc_west ) THEN
757         istop = 0
758         IF( jpiwob < 2 .OR.  jpiwob >= jpiglo ) THEN
759            IF(lwp) WRITE(numout,cform_err)
760            IF(lwp) WRITE(numout,*) '            jpiwob exceed ', jpim1, ' or less than 2'
761            istop = istop + 1
762         END IF
763
764         IF( lk_mpp ) THEN
765            IF( (nimpp < jpiwob+5) .AND. (nimpp > 1) ) THEN
766               IF(lwp) WRITE(numout,cform_err)
767               IF(lwp) WRITE(numout,*) '        A sub-domain is too close to the West OBC'
768               IF(lwp) WRITE(numout,*) '        nimpp must be > jpiwob-5 or =1'
769               istop = istop + 1
770            ENDIF
771         ELSE
772   
773            ! ... stop if  e r r o r (s)   detected
774            IF( istop /= 0 ) THEN
775               WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop'
776               CALL ctl_stop( ctmp1 )
777            ENDIF
778         ENDIF
779      ENDIF
780
781      ! control of the north boundary
782      IF( lp_obc_north ) THEN
783         istop = 0
784         IF( jpjnob < 4 .OR.  jpjnob >= jpjglo ) THEN
785            IF(lwp) WRITE(numout,cform_err)
786            IF(lwp) WRITE(numout,*) '          jpjnob exceed ', jpjm1,' or less than 4'
787            istop = istop + 1
788         END IF
789
790         IF( lk_mpp ) THEN
791            IF( njmpp > jpjnob-5) THEN
792               IF(lwp) WRITE(numout,cform_err)
793               IF(lwp) WRITE(numout,*) '        A sub-domain is too close to the North OBC'
794               IF(lwp) WRITE(numout,*) '        njmpp must be < jpjnob-5'
795               istop = istop + 1
796            ENDIF
797         ELSE
798   
799            ! ... stop if  e r r o r (s)   detected
800            IF( istop /= 0 ) THEN
801                WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop'
802               CALL ctl_stop( ctmp1 )
803           ENDIF
804         ENDIF
805      ENDIF
806
807      ! control of the south boundary
808      IF( lp_obc_south ) THEN
809         istop = 0
810         IF( jpjsob < 2 .OR. jpjsob >= jpjglo ) THEN
811            IF(lwp) WRITE(numout,cform_err)
812            IF(lwp) WRITE(numout,*) '          jpjsob exceed ', jpjm1,' or less than 2'
813            istop = istop + 1
814         END IF
815
816         IF( lk_mpp ) THEN
817            IF( (njmpp < jpjsob+5) .AND. (njmpp > 1) ) THEN
818               IF(lwp) WRITE(numout,cform_err)
819               IF(lwp) WRITE(numout,*) '        A sub-domain is too close to the South OBC'
820               IF(lwp) WRITE(numout,*) '        njmpp must be > jpjsob+5 or =1'
821               istop = istop + 1
822            ENDIF
823         ELSE
824   
825            ! ... stop if  e r r o r (s)   detected
826            IF( istop /= 0 ) THEN
827               WRITE(ctmp1,*) istop,' obcini : E R R O R (S) detected : stop'
828               CALL ctl_stop( ctmp1 )
829            ENDIF
830         ENDIF
831      ENDIF
832
833   END SUBROUTINE obc_init
834
835#else
836   !!---------------------------------------------------------------------------------
837   !!   Dummy module                                                NO open boundaries
838   !!---------------------------------------------------------------------------------
839CONTAINS
840   SUBROUTINE obc_init      ! Dummy routine
841   END SUBROUTINE obc_init
842#endif
843
844   !!=================================================================================
845END MODULE obcini
Note: See TracBrowser for help on using the repository browser.