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

source: trunk/NEMO/OPA_SRC/OBC/obcdta.F90 @ 526

Last change on this file since 526 was 473, checked in by opalod, 18 years ago

nemo_v1_update_060: SM: IOM + 301 levels + CORE + begining of ctl_stop
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 51.8 KB
Line 
1MODULE obcdta
2   !!==============================================================================
3   !!                            ***  MODULE obcdta  ***
4   !! Open boundary data : read the data for the open boundaries.
5   !!==============================================================================
6   !! History :
7   !!        !  98-05 (J.M. Molines) Original code
8   !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90
9   !!   9.0  !  04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input
10   !!------------------------------------------------------------------------------
11#if defined key_obc
12   !!------------------------------------------------------------------------------
13   !!   'key_obc'         :                                Open Boundary Conditions
14   !!------------------------------------------------------------------------------
15   !!   obc_dta           : read u, v, t, s data along each open boundary
16   !!   obc_dta_psi       : read psi data along each open boundary (rigid lid only)
17   !!------------------------------------------------------------------------------
18   !! * Modules used
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
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 daymod          ! calendar
25   USE in_out_manager  ! I/O logical units
26   USE lib_mpp         ! distributed memory computing
27   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
28   USE iom
29#  if defined key_dynspg_rl
30   USE obccli          ! climatological obc, use only in rigid-lid case
31#  endif
32
33   IMPLICIT NONE
34   PRIVATE
35
36   !! * Accessibility
37   PUBLIC obc_dta        ! routines called by step.F90
38   PUBLIC obc_dta_bt     ! routines called by dynspg_ts.F90
39
40   !! * Shared module variables
41   INTEGER ::   &
42      nlecto,   &  ! switch for the first read
43      ntobc1,   &  ! first record used
44      ntobc2,   &  ! second record used
45      ntobc        ! number of time steps in OBC files
46
47   REAL(wp), DIMENSION(:), ALLOCATABLE :: tcobc      ! time_counter variable of BCs
48
49   !! * Substitutions
50#  include "domzgr_substitute.h90"                                           
51#  include "obc_vectopt_loop_substitute.h90"
52   !!---------------------------------------------------------------------------------
53   !!   OPA 9.0 , LODYC-IPSL  (2003)
54   !! $Header$
55   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
56   !!---------------------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE obc_dta( kt )
61      !!--------------------------------------------------------------------
62      !!              ***  SUBROUTINE obc_dta  ***
63      !!                   
64      !! ** Purpose :   Find the climatological boundary arrays for the specified date,
65      !!   The boundary arrays are netcdf files. Three possible cases:
66      !!   - one time frame only in the file (time dimension = 1).
67      !!     in that case the boundary data does not change in time.
68      !!   - many time frames. In that case,  if we have 12 frames
69      !!     we assume monthly fields.
70      !!     Else, we assume that time_counter is in seconds
71      !!     since the beginning of either the current year or a reference
72      !!     year given in the namelist.
73      !!     (no check is done so far but one would have to check the "unit"
74      !!     attribute of variable time_counter).
75      !!
76      !! History :
77      !!        !  98-05 (J.M. Molines) Original code
78      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90
79      !!   9.0  !  04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input
80      !!--------------------------------------------------------------------
81      !! * Arguments
82      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
83     
84      !! * Local declarations
85      INTEGER ::   ji, jj, jk, ii, ij   ! dummy loop indices
86      INTEGER ::   itimo, iman, imois
87      INTEGER ::   i15
88      REAL(wp) ::   zxy
89      !! * Ajouts FD
90      INTEGER ::  isrel              ! number of seconds since 1/1/1992
91      INTEGER, DIMENSION(1) ::  itobce, itobcw,  & ! number of time steps in OBC files
92                                itobcs, itobcn     !    "       "       "       "
93      INTEGER ::  istop       
94      INTEGER ::  iprint        ! frequency for printouts.
95      INTEGER :: id_e, id_w, id_n, id_s       ! file identifiers
96      LOGICAL :: llnot_done
97      CHARACTER(LEN=25) :: cl_vname
98      !!--------------------------------------------------------------------
99
100      IF( lk_dynspg_rl )  THEN
101         CALL obc_dta_psi (kt)     ! update bsf data at open boundaries
102         IF ( nobc_dta == 1 .AND. kt == nit000 ) CALL ctl_stop( 'obcdta : time-variable psi boundary data not allowed yet' )
103      ENDIF
104           
105      ! 1.   First call: check time frames available in files.
106      ! -------------------------------------------------------
107
108      IF ( kt == nit000 ) THEN
109     
110         nlecto =  0
111
112         IF (lwp) WRITE(numout,*)
113         IF (lwp) WRITE(numout,*)     'obc_dta : find boundary data'
114         IF (lwp) WRITE(numout,*)     '~~~~~~~'
115             
116         IF ( nobc_dta == 0 ) THEN
117            IF(lwp) WRITE(numout,*)  '  OBC data taken from initial conditions.'
118            ntobc1 = 1
119            ntobc2 = 1
120         ELSE   
121            IF (lwp) WRITE(numout,*)  '  OBC data taken from netcdf files.'
122            IF (lwp) WRITE(numout,*)  '  climatology (T/F):',ln_obc_clim
123            ! check the number of time steps in the files.
124            cl_vname = 'time_counter'
125            IF ( lp_obc_east ) THEN
126               CALL iom_open ( 'obceast_TS.nc' , id_e )
127               idvar = iom_varid( id_e, cl_vname, kdimsz = itobce )
128            ENDIF
129            IF ( lp_obc_west ) THEN
130               CALL iom_open ( 'obcwest_TS.nc' , id_w )
131               idvar = iom_varid( id_w, cl_vname, kdimsz = itobcw )
132            ENDIF
133            IF ( lp_obc_north ) THEN
134               CALL iom_open ( 'obcnorth_TS.nc', id_n )
135               idvar = iom_varid( id_n, cl_vname, kdimsz = itobcn )
136            ENDIF
137            IF ( lp_obc_south ) THEN
138               CALL iom_open ( 'obcsouth_TS.nc', id_s )
139               idvar = iom_varid( id_s, cl_vname, kdimsz = itobcs )
140            ENDIF
141
142            ntobc = MAX(itobce(1),itobcw(1),itobcn(1),itobcs(1))
143            istop = 0
144            IF ( lp_obc_east  .AND. itobce(1) /= ntobc ) istop = 1 
145            IF ( lp_obc_west  .AND. itobcw(1) /= ntobc ) istop = 1     
146            IF ( lp_obc_north .AND. itobcn(1) /= ntobc ) istop = 1
147            IF ( lp_obc_south .AND. itobcs(1) /= ntobc ) istop = 1 
148            IF ( istop /= 0 )  THEN
149               WRITE(ctmp1,*) ' east, west, north, south: ', itobce(1), itobcw(1), itobcn(1), itobcs(1)
150               CALL ctl_stop( 'obcdta : all files must have the same number of time steps', ctmp1 )
151            ENDIF
152            IF ( ntobc == 1 ) THEN
153               IF ( lwp ) WRITE(numout,*) ' obcdta found one time step only in the OBC files'
154            ELSE
155               ALLOCATE (tcobc(ntobc))
156               llnot_done = .TRUE.
157               IF ( lp_obc_east ) THEN
158                  IF ( llnot_done ) THEN
159                     CALL iom_gettime ( id_e, TRIM(cl_vname), tcobc )
160                     llnot_done = .FALSE.
161                  ENDIF
162                  CALL iom_close (id_e)
163               ENDIF
164               IF ( lp_obc_west ) THEN
165                  IF ( llnot_done ) THEN
166                     CALL iom_gettime ( id_w, TRIM(cl_vname), tcobc )
167                     llnot_done = .FALSE.
168                 ENDIF
169                 CALL iom_close (id_w)
170               ENDIF
171               IF ( lp_obc_north ) THEN
172                  IF ( llnot_done ) THEN
173                     CALL iom_gettime ( id_n, TRIM(cl_vname), tcobc )
174                     llnot_done = .FALSE.
175                  ENDIF
176                  CALL iom_close (id_n)
177               ENDIF
178               IF ( lp_obc_south ) THEN
179                  IF ( llnot_done ) THEN
180                     CALL iom_gettime ( id_s, TRIM(cl_vname), tcobc )
181                     llnot_done = .FALSE.
182                  ENDIF
183                  CALL iom_close (id_s)
184               ENDIF
185               IF ( lwp ) WRITE(numout,*) ' obcdta found', ntobc,' time steps in the OBC files'
186               IF ( .NOT. ln_obc_clim .AND. ntobc == 12 ) THEN
187                  IF ( lwp ) WRITE(numout,*) '  WARNING: With monthly data we assume climatology'
188                  ln_obc_clim = .true.
189               ENDIF
190            ENDIF
191         ENDIF
192     
193       ! 1.1  Tangential velocities set to zero
194       ! --------------------------------------
195         IF( lp_obc_east  ) vfoe = 0.e0
196         IF( lp_obc_west  ) vfow = 0.e0
197         IF( lp_obc_south ) ufos = 0.e0
198         IF( lp_obc_north ) ufon = 0.e0
199     
200       ! 1.2  Data temperature, salinity, normal velocities set to zero
201       !                        or initial conditions if nobc_dta == 0
202       ! --------------------------------------------------------------
203
204         IF( lp_obc_east )   THEN
205            ! initialisation to zero
206            sedta(:,:,:) = 0.e0
207            tedta(:,:,:) = 0.e0
208            uedta(:,:,:) = 0.e0
209            !                                    ! ================== !
210            IF( nobc_dta == 0 )   THEN           ! initial state used
211               !                                 ! ================== !
212               !  Fills sedta, tedta, uedta (global arrays)
213               !  Remark: this works for njzoom = 1.
214               !          Should the definition of ij include njzoom?
215               DO ji = nie0, nie1
216                  DO jk = 1, jpkm1
217                     DO jj = nje0p1, nje1m1
218                        ij = jj -1 + njmpp
219                        sedta(ij,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk)
220                        tedta(ij,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk)
221                        uedta(ij,jk,1) = un(ji,jj,jk)*umask(ji,jj,jk)
222                     END DO
223                  END DO
224               END DO
225            ENDIF
226         ENDIF
227
228         IF( lp_obc_west )   THEN
229            ! initialisation to zero
230            swdta(:,:,:) = 0.e0
231            twdta(:,:,:) = 0.e0
232            uwdta(:,:,:) = 0.e0
233            !                                    ! ================== !
234            IF( nobc_dta == 0 )   THEN           ! initial state used !
235               !                                 ! ================== !
236               !  Fills swdta, twdta, uwdta (global arrays)
237               !  Remark: this works for njzoom = 1.
238               !          Should the definition of ij include njzoom?
239               DO ji = niw0, niw1
240                  DO jk = 1, jpkm1
241                     DO jj = njw0p1, njw1m1
242                        ij = jj -1 + njmpp
243                        swdta(ij,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk)
244                        twdta(ij,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk)
245                        uwdta(ij,jk,1) = un(ji,jj,jk)*umask(ji,jj,jk)
246                     END DO
247                  END DO
248               END DO
249            ENDIF
250         ENDIF
251
252         IF( lp_obc_north)   THEN
253            ! initialisation to zero
254            sndta(:,:,:) = 0.e0
255            tndta(:,:,:) = 0.e0
256            vndta(:,:,:) = 0.e0
257            !                                    ! ================== !
258            IF( nobc_dta == 0 )   THEN           ! initial state used
259               !                                 ! ================== !
260               !  Fills sndta, tndta, vndta (global arrays)
261               !  Remark: this works for njzoom = 1.
262               !          Should the definition of ij include njzoom?
263               DO jj = njn0, njn1
264                  DO jk = 1, jpkm1
265                     DO ji = nin0p1, nin1m1
266                        ii = ji -1 + nimpp
267                        sndta(ii,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk)
268                        tndta(ii,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk)
269                        vndta(ii,jk,1) = vn(ji,jj,jk)*vmask(ji,jj,jk)
270                     END DO
271                  END DO
272               END DO
273            ENDIF
274         ENDIF
275
276         IF( lp_obc_south )   THEN
277            ! initialisation to zero
278            ssdta(:,:,:) = 0.e0
279            tsdta(:,:,:) = 0.e0
280            vsdta(:,:,:) = 0.e0
281            !                                    ! ================== !
282            IF( nobc_dta == 0 )   THEN           ! initial state used
283               !                                 ! ================== !
284               !  Fills ssdta, tsdta, vsdta (global arrays)
285               !  Remark: this works for njzoom = 1.
286               !          Should the definition of ij include njzoom?
287               DO jj = njs0, njs1
288                  DO jk = 1, jpkm1
289                     DO ji = nis0p1, nis1m1
290                        ii = ji -1 + nimpp
291                        ssdta(ii,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk)
292                        tsdta(ii,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk)
293                        vsdta(ii,jk,1) = vn(ji,jj,jk)*vmask(ji,jj,jk)
294                     END DO
295                  END DO
296               END DO
297            ENDIF
298         ENDIF
299
300      ENDIF        !       end if kt == nit000
301     
302      ! 2.  Initialize the time we are at.
303      !     Does this every time the routine is called,
304      !     excepted when nobc_dta = 0
305      !---------------------------------------------------------------------
306      IF( nobc_dta == 0 )   THEN
307         itimo = 1
308         zxy   = 0
309      ELSE
310         IF( ntobc == 1 )   THEN
311            itimo = 1
312         ELSE IF( ntobc == 12 )   THEN      !   BC are monthly   
313            ! we assume we have climatology in that case
314            iman  = 12
315            i15   = nday / 16
316            imois = nmonth + i15 - 1
317            IF( imois == 0 )   imois = iman
318            itimo = imois   
319         ELSE
320            IF(lwp) WRITE(numout,*) 'data other than constant or monthly', kt
321            iman  = ntobc
322            itimo = FLOOR( kt*rdt / (tcobc(2)-tcobc(1)) )
323            isrel = kt*rdt
324         ENDIF
325      ENDIF
326     
327      ! 2.1 Read two records in the file if necessary
328      ! ---------------------------------------------
329      IF( ( nobc_dta == 1 ) .AND. ( ( kt == nit000 .AND. nlecto == 0 ) .OR. itimo  /= ntobc1 ) )   THEN
330         nlecto = 1
331     
332         ! Calendar computation
333         IF( ntobc == 1 )   THEN            !  BC are constant in time
334            ntobc1 = 1
335            ntobc2 = 1 
336         ELSE IF( ntobc == 12 )   THEN      !   BC are monthly   
337            ntobc1 = itimo         ! first file record used
338            ntobc2 = ntobc1 + 1    ! last  file record used
339            ntobc1 = MOD( ntobc1, iman )
340            IF( ntobc1 == 0 )   ntobc1 = iman
341            ntobc2 = MOD( ntobc2, iman )
342            IF( ntobc2 == 0 )   ntobc2 = iman
343            IF( lwp )   THEN
344               WRITE(numout,*) ' read monthly obc first record file used ntobc1 ', ntobc1
345               WRITE(numout,*) ' read monthly obc last  record file used ntobc2 ', ntobc2
346            ENDIF
347         ELSE
348            isrel=kt*rdt
349            ntobc1 = itimo         ! first file record used
350            ntobc2 = ntobc1 + 1    ! last  file record used
351            ntobc1 = MOD( ntobc1, iman )
352            IF( ntobc1 == 0 )   ntobc1 = iman
353            ntobc2 = MOD( ntobc2, iman )
354            IF( ntobc2 == 0 )   ntobc2 = iman
355            IF(lwp) WRITE(numout,*) ' read obc first record file used ntobc1 ', ntobc1
356            IF(lwp) WRITE(numout,*) ' read obc last  record file used ntobc2 ', ntobc2
357         ENDIF
358                               ! ======================= !
359                               !  BCs read               !
360         !                     ! ======================= !
361         IF( lp_obc_east )   THEN
362            ! ... Read datafile and set temperature, salinity and normal velocity
363            ! ... initialise the sedta, tedta, uedta arrays
364            CALL iom_open ( 'obceast_TS.nc' , id_e )
365            CALL iom_get ( id_e, jpdom_unknown, 'votemper', tedta(:,:,1), ktime=ntobc1 )
366            CALL iom_get ( id_e, jpdom_unknown, 'votemper', tedta(:,:,2), ktime=ntobc2 )
367            CALL iom_get ( id_e, jpdom_unknown, 'vosaline', sedta(:,:,1), ktime=ntobc1 )
368            CALL iom_get ( id_e, jpdom_unknown, 'vosaline', sedta(:,:,2), ktime=ntobc2 )
369            CALL iom_close (id_e)
370            !
371            CALL iom_open ( 'obceast_U.nc' , id_e )
372            CALL iom_get ( id_e, jpdom_unknown, 'vozocrtx', uedta(:,:,1), ktime=ntobc1 )
373            CALL iom_get ( id_e, jpdom_unknown, 'vozocrtx', uedta(:,:,2), ktime=ntobc2 )
374            CALL iom_close ( id_e )
375            !  Usually printout is done only once at kt = nit000,
376            !  unless nprint (namelist) > 1     
377            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN
378               WRITE(numout,*)
379               WRITE(numout,*) ' Read East OBC data records ', ntobc1, ntobc2
380               iprint = (jpjef-jpjed+1)/20 +1
381               WRITE(numout,*) ' Temperature record 1 - printout every 3 level'
382               CALL prihre( tedta(:,:,1),jpjef-jpjed+1,jpk,1,jpjef-jpjed+1,iprint, &
383                  &         jpk, 1, -3, 1., numout )
384               WRITE(numout,*)
385               WRITE(numout,*) ' Salinity  record 1 - printout every 3 level'
386               CALL prihre( sedta(:,:,1), jpjef-jpjed+1, jpk, 1, jpjef-jpjed+1, iprint, &
387                  &        jpk, 1, -3, 1., numout )
388               WRITE(numout,*)
389               WRITE(numout,*) ' Normal velocity U  record 1  - printout every 3 level'
390               CALL prihre( uedta(:,:,1), jpjef-jpjed+1, jpk, 1, jpjef-jpjed+1, iprint, &
391                  &         jpk, 1, -3, 1., numout )
392            ENDIF
393         ENDIF
394
395         IF( lp_obc_west )   THEN
396            ! ... Read datafile and set temperature, salinity and normal velocity
397            ! ... initialise the swdta, twdta, uwdta arrays
398            CALL iom_open ( 'obcwest_TS.nc' , id_w )
399            CALL iom_get ( id_w, jpdom_unknown, 'votemper', twdta(:,:,1), ktime=ntobc1 )
400            CALL iom_get ( id_w, jpdom_unknown, 'votemper', twdta(:,:,2), ktime=ntobc2 )
401            CALL iom_get ( id_w, jpdom_unknown, 'vosaline', swdta(:,:,1), ktime=ntobc1 )
402            CALL iom_get ( id_w, jpdom_unknown, 'vosaline', swdta(:,:,2), ktime=ntobc2 )
403            CALL iom_close (id_w)
404            !
405            CALL iom_open ( 'obcwest_U.nc' , id_w )
406            CALL iom_get ( id_w, jpdom_unknown, 'vozocrtx', uwdta(:,:,1), ktime=ntobc1 )
407            CALL iom_get ( id_w, jpdom_unknown, 'vozocrtx', uwdta(:,:,2), ktime=ntobc2 )
408            CALL iom_close ( id_w )
409            !
410            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN
411               WRITE(numout,*)
412               WRITE(numout,*) ' Read West OBC data records ', ntobc1, ntobc2
413               iprint = (jpjwf-jpjwd+1)/20 +1
414               WRITE(numout,*) ' Temperature  record 1 - printout every 3 level'
415               CALL prihre( twdta(:,:,1), jpjwf-jpjwd+1, jpk, 1, jpjwf-jpjwd+1, iprint, &
416                  &         jpk, 1, -3, 1., numout )
417               WRITE(numout,*)
418               WRITE(numout,*) ' Salinity  record 1 - printout every 3 level'
419               CALL prihre( swdta(:,:,1), jpjwf-jpjwd+1, jpk, 1, jpjwf-jpjwd+1, iprint, &
420                  &         jpk, 1, -3, 1., numout )
421               WRITE(numout,*)
422               WRITE(numout,*) ' Normal velocity U  record 1  - printout every 3 level'
423               CALL prihre( uwdta(:,:,1), jpjwf-jpjwd+1, jpk, 1, jpjwf-jpjwd+1, iprint, &
424                  &         jpk, 1, -3, 1., numout )
425            ENDIF
426         ENDIF
427
428         IF( lp_obc_north )   THEN     
429            CALL iom_open ( 'obcnorth_TS.nc', id_n )
430            CALL iom_get ( id_n, jpdom_unknown, 'votemper', tndta(:,:,1), ktime=ntobc1 )
431            CALL iom_get ( id_n, jpdom_unknown, 'votemper', tndta(:,:,2), ktime=ntobc2 )
432            CALL iom_get ( id_n, jpdom_unknown, 'vosaline', sndta(:,:,1), ktime=ntobc1 )
433            CALL iom_get ( id_n, jpdom_unknown, 'vosaline', sndta(:,:,2), ktime=ntobc2 )
434            CALL iom_close ( id_n )                                                           
435            !
436            CALL iom_open ( 'obcnorth_V.nc', id_n )                                         
437            CALL iom_get ( id_n, jpdom_unknown, 'vomecrty', vndta(:,:,1), ktime=ntobc1 )
438            CALL iom_get ( id_n, jpdom_unknown ,'vomecrty', vndta(:,:,2), ktime=ntobc2 )
439            CALL iom_close ( id_n )
440            !
441            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN
442               WRITE(numout,*)
443               WRITE(numout,*) ' Read North OBC data records ', ntobc1, ntobc2
444               iprint = (jpinf-jpind+1)/20 +1
445               WRITE(numout,*) ' Temperature  record 1 - printout every 3 level'
446               CALL prihre( tndta(:,:,1), jpinf-jpind+1, jpk, 1, jpinf-jpind+1, iprint, &
447                  &         jpk, 1, -3, 1., numout )
448               WRITE(numout,*)
449               WRITE(numout,*) ' Salinity  record 1 - printout every 3 level'
450               CALL prihre( sndta(:,:,1), jpinf-jpind+1, jpk, 1, jpinf-jpind+1, iprint, &
451                  &         jpk, 1, -3, 1., numout )
452               WRITE(numout,*)
453               WRITE(numout,*) ' Normal velocity V  record 1  - printout every 3 level'
454               CALL prihre( vndta(:,:,1), jpinf-jpind+1, jpk, 1, jpinf-jpind+1, iprint, &
455                  &         jpk, 1, -3, 1., numout )
456            ENDIF
457         ENDIF
458
459         IF( lp_obc_south )   THEN     
460            CALL iom_open ( 'obcsouth_TS.nc', id_s )
461            CALL iom_get ( id_s, jpdom_unknown, 'votemper', tsdta(:,:,1), ktime=ntobc1 )
462            CALL iom_get ( id_s, jpdom_unknown, 'votemper', tsdta(:,:,2), ktime=ntobc2 )
463            CALL iom_get ( id_s, jpdom_unknown, 'vosaline', ssdta(:,:,1), ktime=ntobc1 )
464            CALL iom_get ( id_s, jpdom_unknown, 'vosaline', ssdta(:,:,2), ktime=ntobc2 )
465            CALL iom_close ( id_s )                                                           
466            !
467            CALL iom_open ( 'obcsouth_V.nc', id_s )                                         
468            CALL iom_get ( id_s, jpdom_unknown, 'vomecrty', vsdta(:,:,1), ktime=ntobc1 )
469            CALL iom_get ( id_s, jpdom_unknown ,'vomecrty', vsdta(:,:,2), ktime=ntobc2 )
470            CALL iom_close ( id_s )
471            !
472            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 ) )   THEN
473               WRITE(numout,*)
474               WRITE(numout,*) ' Read South OBC data records ', ntobc1, ntobc2
475               iprint = (jpisf-jpisd+1)/20 +1
476               WRITE(numout,*) ' Temperature  record 1 - printout every 3 level'
477               CALL prihre( tsdta(:,:,1), jpisf-jpisd+1, jpk, 1, jpisf-jpisd+1, iprint, &
478                  &         jpk, 1, -3, 1., numout )
479               WRITE(numout,*)
480               WRITE(numout,*) ' Salinity  record 1 - printout every 3 level'
481               CALL prihre( ssdta(:,:,1), jpisf-jpisd+1, jpk, 1, jpisf-jpisd+1, iprint, &
482                  &         jpk, 1, -3, 1., numout )
483               WRITE(numout,*)
484               WRITE(numout,*) ' Normal velocity V  record 1  - printout every 3 level'
485               CALL prihre( vsdta(:,:,1), jpisf-jpisd+1, jpk, 1, jpisf-jpisd+1, iprint, &
486                  &         jpk, 1, -3, 1., numout )
487            ENDIF
488         ENDIF
489
490      ELSE
491         
492         nlecto = 0        !      no reading of OBC barotropic data                         
493
494      ENDIF                !      end of the test on the condition to read or not the files
495     
496      ! 3.  Call at every time step :
497      !     Linear interpolation of BCs to current time step
498      ! ----------------------------------------------------
499
500      IF( ntobc == 1 .OR. nobc_dta == 0 )   THEN
501         zxy = 0.
502      ELSE IF( ntobc == 12 )   THEN         
503         zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
504      ELSE
505         zxy = (tcobc(ntobc1)-FLOAT(isrel))/(tcobc(ntobc1)-tcobc(ntobc2))
506      ENDIF
507     
508      IF( lp_obc_east )   THEN
509         !  fills sfoe, tfoe, ufoe (local to each processor)
510         DO jk = 1, jpkm1
511            DO jj = nje0p1, nje1m1
512               ij = jj -1 + njmpp
513               sfoe(jj,jk) =  ( zxy * sedta(ij,jk,2) + &
514                  &           (1.-zxy) * sedta(ij,jk,1) ) * temsk(jj,jk)
515               tfoe(jj,jk) =  ( zxy * tedta(ij,jk,2) + &
516                  &           (1.-zxy) * tedta(ij,jk,1) ) * temsk(jj,jk)
517               ufoe(jj,jk) =  ( zxy * uedta(ij,jk,2) + &
518                  &           (1.-zxy) * uedta(ij,jk,1) ) * uemsk(jj,jk)
519            END DO
520         END DO
521      ENDIF
522     
523      IF( lp_obc_west )   THEN
524         !  fills sfow, tfow, ufow (local to each processor)
525         DO jk = 1, jpkm1
526            DO jj = njw0p1, njw1m1
527               ij = jj -1 + njmpp
528               sfow(jj,jk) =  ( zxy * swdta(ij,jk,2) + &
529                  &           (1.-zxy) * swdta(ij,jk,1) ) * twmsk(jj,jk)
530               tfow(jj,jk) =  ( zxy * twdta(ij,jk,2) + &
531                  &           (1.-zxy) * twdta(ij,jk,1) ) * twmsk(jj,jk)
532               ufow(jj,jk) =  ( zxy * uwdta(ij,jk,2) + &
533                  &           (1.-zxy) * uwdta(ij,jk,1) ) * uwmsk(jj,jk)
534            END DO
535         END DO
536      ENDIF
537     
538      IF( lp_obc_north )   THEN
539         !  fills sfon, tfon, vfon (local to each processor)
540         DO jk = 1, jpkm1
541            DO ji = nin0p1, nin1m1
542               ii = ji -1 + nimpp
543               sfon(ji,jk) =  ( zxy * sndta(ii,jk,2) + &
544                  &           (1.-zxy) * sndta(ii,jk,1) ) * tnmsk(ji,jk)
545               tfon(ji,jk) =  ( zxy * tndta(ii,jk,2) + &
546                  &           (1.-zxy) * tndta(ii,jk,1) ) * tnmsk(ji,jk)
547               vfon(ji,jk) =  ( zxy * vndta(ii,jk,2) + &
548                  &           (1.-zxy) * vndta(ii,jk,1) ) * vnmsk(ji,jk)
549            END DO
550         END DO
551      ENDIF
552     
553      IF( lp_obc_south )   THEN
554         !  fills sfos, tfos, vfos (local to each processor)
555         DO jk = 1, jpkm1
556           DO ji = nis0p1, nis1m1
557              ii = ji -1 + nimpp
558              sfos(ji,jk) = ( zxy * ssdta(ii,jk,2) + &
559                 &          (1.-zxy) * ssdta(ii,jk,1) ) * tsmsk(ji,jk)
560              tfos(ji,jk) = ( zxy * tsdta(ii,jk,2) + &
561                 &          (1.-zxy) * tsdta(ii,jk,1) ) * tsmsk(ji,jk)
562              vfos(ji,jk) = ( zxy * vsdta(ii,jk,2) + &
563                 &          (1.-zxy) * vsdta(ii,jk,1) ) * vsmsk(ji,jk)
564           END DO
565         END DO             
566      ENDIF
567     
568   END SUBROUTINE obc_dta
569     
570# if defined key_dynspg_rl
571   !!-----------------------------------------------------------------------------
572   !!   Rigid-lid
573   !!-----------------------------------------------------------------------------
574
575   SUBROUTINE obc_dta_psi ( kt )
576      !!-----------------------------------------------------------------------------
577      !!                       ***  SUBROUTINE obc_dta_psi  ***
578      !!
579      !! ** Purpose :
580      !!      Update the climatological streamfunction OBC at each time step.
581      !!      Depends on the user's configuration.  Here data are read only once
582      !!      at the beginning of the run.
583      !!
584      !! ** Method :
585      !!      1. initialization
586      !!         kbsfstart: number of time steps over which increase bsf
587      !!         during initialization. This is provided for a smooth start
588      !!         in cases where the transport is large (like on the Antarctic
589      !!         continent). also note that when kbfstart=1, the transport
590      !!         increases a lot in one time step and the precision usually
591      !!         required for the solver may not be enough.
592      !!      2. set the time evolution of the climatological barotropic streamfunction
593      !!         along the isolated coastlines ( gcbic(jnic) ).
594      !!      3. set the climatological barotropic streamfunction at the boundary.
595      !!
596      !!      The last two steps are done only at first step (nit000) or if kt <= kbfstart
597      !!
598      !! History :
599      !!        ! 97-08 (G. Madec, J.M. Molines)
600      !!   8.5  ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
601      !!   9.0  ! 05-11  (V. Garnier) Surface pressure gradient organization
602      !!----------------------------------------------------------------------------
603      !! * Arguments
604      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
605
606      !! * Local declarations
607      INTEGER ::   ji, jj, jnic, jip         ! dummy loop indices
608      INTEGER ::   inum = 11                 ! temporary logical unit
609      INTEGER ::   ip, ii, ij, iii, ijj
610      INTEGER ::   kbsfstart
611      REAL(wp) ::   zsver1, zsver2, zsver3, z2dtr, zcoef
612      !!----------------------------------------------------------------------------
613
614      ! 1. initialisation
615      ! -----------------
616
617      kbsfstart =  1
618      zsver1 =  bsfic0(1)
619      zsver2 =  zsver1
620      IF( kt <= kbsfstart )   THEN
621         zcoef = float(kt)/float(kbsfstart)
622      ELSE
623         zcoef = 1.
624      END IF
625      bsfic(1) = zsver1*zcoef
626      IF( lwp .AND. ( kt <= kbsfstart ) )   THEN
627         IF(lwp) WRITE(numout,*)'            '
628         IF(lwp) WRITE(numout,*)'obcdta: spinup phase in obc_dta_psi routine'
629         IF(lwp) WRITE(numout,*)'~~~~~~  it=',kt,'  OBC: spinup coef: ', &
630                                           zcoef, ' and transport: ',bsfic(1)
631      END IF
632
633      zsver2 =  bsfic(1)-bsfic(2)
634      zsver3 =  bsfic(2)
635
636      ! 2. Right hand side of the barotropic elliptic equation (isolated coastlines)
637      ! ----------------------------------------------------------------------------
638
639      IF( ( neuler == 0 ) .AND. ( kt == nit000 ) )   THEN
640         z2dtr = 1./rdt
641      ELSE
642         z2dtr = 1./2./rdt
643      END IF
644      ! ... bsfb(ii,ij) should be constant but due to the Asselin filter it
645      ! ... converges asymptotically towards bsfic(jnic)
646      ! ... However, bsfb(ii,ij) is constant along the same coastlines
647      ! ... ---> can be improved using an extra array for storing bsficb (before)
648      IF( nbobc > 1 )   THEN
649         DO jnic = 1,nbobc - 1
650            gcbic(jnic) = 0.e0
651            ip=mnic(0,jnic)
652            DO jip = 1,ip
653               ii = miic(jip,0,jnic) 
654               ij = mjic(jip,0,jnic) 
655               IF( ii >= nldi+ nimpp - 1 .AND. ii <= nlei+ nimpp - 1 .AND. &
656                   ij >= nldj+ njmpp - 1 .AND. ij <= nlej+ njmpp - 1 )   THEN
657                  iii=ii-nimpp+1
658                  ijj=ij-njmpp+1
659                  gcbic(jnic) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
660               END IF
661            END DO
662         END DO
663      END IF
664
665      IF( lk_mpp )   CALL mpp_isl( gcbic, 3 )
666
667      ! 3. Update the climatological barotropic function at the boundary
668      ! ----------------------------------------------------------------
669
670      IF( lpeastobc )   THEN
671
672         IF( kt == nit000 .OR. kt <= kbsfstart )   THEN
673            OPEN(inum,file='obceastbsf.dta')
674            READ(inum,*)
675            READ(inum,*)
676            READ(inum,*)
677            READ(inum,*)
678            READ(inum,*)
679            READ(inum,*) (bfoe(jj),jj=jpjed, jpjef)
680            CLOSE(inum)
681         END IF
682         DO jj=jpjed, jpjefm1
683            bfoe(jj)=bfoe(jj)*zcoef
684         END DO
685
686      END IF
687
688      IF( lpwestobc)   THEN
689
690         IF( kt == nit000 .OR. kt <= kbsfstart ) then
691            OPEN(inum,file='obcwestbsf.dta')
692            READ(inum,*)
693            READ(inum,*)
694            READ(inum,*)
695            READ(inum,*)
696            READ(inum,*)
697            READ(inum,*) (bfow(jj),jj=jpjwd, jpjwf)
698            CLOSE(inum)
699         END IF
700         DO jj=jpjwd, jpjwfm1
701            bfow(jj)=bfow(jj)*zcoef
702         END DO
703
704      END IF
705
706      IF( lpsouthobc)   THEN
707
708         IF( kt == nit000 .OR. kt <= kbsfstart )   THEN
709            OPEN(inum,file='obcsouthbsf.dta')
710            READ(inum,*)
711            READ(inum,*)
712            READ(inum,*)
713            READ(inum,*)
714            READ(inum,*)
715            READ(inum,*) (bfos(jj),jj=jpisd, jpisf)
716            CLOSE(inum)
717         END IF
718         DO ji=jpisd, jpisfm1
719            bfos(ji)=bfos(ji)*zcoef
720         END DO
721
722      END IF
723
724      IF( lpnorthobc)   THEN
725         IF( kt == nit000 .OR. kt <= kbsfstart )   THEN
726            OPEN(inum,file='obcnorthbsf.dta')
727            READ(inum,*)
728            READ(inum,*)
729            READ(inum,*)
730            READ(inum,*)
731            READ(inum,*)
732            READ(inum,*) (bfon(jj),jj=jpind, jpinf)
733            CLOSE(inum)
734         END IF
735         DO ji=jpind, jpinfm1
736            bfon(ji)=bfon(ji)*zcoef
737         END DO
738
739      END IF
740
741   END SUBROUTINE obc_dta_psi
742#else
743   !!-----------------------------------------------------------------------------
744   !!   Default option                   
745   !!-----------------------------------------------------------------------------
746   SUBROUTINE obc_dta_psi ( kt )       ! Empty routine
747      !! * Arguments
748      INTEGER,INTENT(in) :: kt 
749      WRITE(*,*) 'obc_dta_psi: You should not have seen this print! error?', kt
750   END SUBROUTINE obc_dta_psi
751# endif
752
753
754#if defined key_dynspg_ts || defined key_dynspg_exp
755   SUBROUTINE obc_dta_bt( kt, kbt )
756      !!---------------------------------------------------------------------------
757      !!                      ***  SUBROUTINE obc_dta  ***
758      !!
759      !! ** Purpose :   time interpolation of barotropic data for time-splitting scheme
760      !!                Data at the boundary must be in m2/s
761      !!
762      !! History :
763      !!   9.0  !  05-11 (V. garnier) Original code
764      !!---------------------------------------------------------------------------
765      !! * Arguments
766      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
767      INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index
768
769      !! * Local declarations
770      INTEGER ::   ji, jj, jk, ii, ij   ! dummy loop indices
771      INTEGER ::   id_e, id_w, id_n, id_s, fid  ! file identifiers
772      INTEGER ::   itimo, iman, imois, i15
773      INTEGER ::   itobcm, itobcp, itimom, itimop
774      REAL(wp) ::  zxy
775      INTEGER ::   isrel, ikt           ! number of seconds since 1/1/1992
776      INTEGER ::   iprint              ! frequency for printouts.
777
778      !!---------------------------------------------------------------------------
779
780      ! 1.   First call: check time frames available in files.
781      ! -------------------------------------------------------
782
783      IF( kt == nit000 ) THEN
784
785         ! 1.1  Barotropic tangential velocities set to zero
786         ! -------------------------------------------------
787         IF( lp_obc_east  ) vbtfoe(:) = 0.e0
788         IF( lp_obc_west  ) vbtfow(:) = 0.e0
789         IF( lp_obc_south ) ubtfos(:) = 0.e0
790         IF( lp_obc_north ) ubtfon(:) = 0.e0
791
792         ! 1.2  Sea surface height and normal barotropic velocities set to zero
793         !                               or initial conditions if nobc_dta == 0
794         ! --------------------------------------------------------------------
795
796          IF( lp_obc_east ) THEN
797             ! initialisation to zero
798             sshedta(:,:) = 0.e0
799             ubtedta(:,:) = 0.e0
800             !                                        ! ================== !
801             IF( nobc_dta == 0 )   THEN               ! initial state used !
802                !                                     ! ================== !
803                !  Fills sedta, tedta, uedta (global arrays)
804                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
805                DO ji = nie0, nie1
806                   DO jj = nje0p1, nje1m1
807                      ij = jj -1 + njmpp
808                      sshedta(ij,1) = sshn(ji+1,jj) * tmask(ji+1,jj,1)
809                   END DO
810                END DO
811             ENDIF
812          ENDIF
813
814          IF( lp_obc_west) THEN
815             ! initialisation to zero
816             sshwdta(:,:) = 0.e0
817             ubtwdta(:,:) = 0.e0
818             !                                        ! ================== !
819             IF( nobc_dta == 0 )   THEN               ! initial state used !
820                !                                     ! ================== !
821                !  Fills swdta, twdta, uwdta (global arrays)
822                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
823                DO ji = niw0, niw1
824                   DO jj = njw0p1, njw1m1
825                      ij = jj -1 + njmpp
826                      sshwdta(ij,1) = sshn(ji,jj) * tmask(ji,jj,1)
827                   END DO
828                END DO
829             ENDIF
830          ENDIF
831
832          IF( lp_obc_north) THEN
833             ! initialisation to zero
834             sshndta(:,:) = 0.e0
835             vbtndta(:,:) = 0.e0
836             !                                        ! ================== !
837             IF( nobc_dta == 0 )   THEN               ! initial state used !
838                !                                     ! ================== !
839                !  Fills sndta, tndta, vndta (global arrays)
840                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
841                DO jj = njn0, njn1
842                   DO ji = nin0p1, nin1m1
843                      DO jk = 1, jpkm1
844                         ii = ji -1 + nimpp
845                         vbtndta(ii,1) = vbtndta(ii,1) + vndta(ii,jk,1)*fse3v(ji,jj,jk)
846                      END DO
847                      sshndta(ii,1) = sshn(ii,jj+1) * tmask(ji,jj+1,1)
848                   END DO
849                END DO
850             ENDIF
851          ENDIF
852
853          IF( lp_obc_south) THEN
854             ! initialisation to zero
855             ssdta(:,:,:) = 0.e0
856             tsdta(:,:,:) = 0.e0
857             vsdta(:,:,:) = 0.e0
858             sshsdta(:,:) = 0.e0
859             vbtsdta(:,:) = 0.e0
860             !                                        ! ================== !
861             IF( nobc_dta == 0 )   THEN               ! initial state used !
862                !                                     ! ================== !
863                !  Fills ssdta, tsdta, vsdta (global arrays)
864                !  Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
865                DO jj = njs0, njs1
866                   DO ji = nis0p1, nis1m1
867                      DO jk = 1, jpkm1
868                         ii = ji -1 + nimpp
869                         vbtsdta(ii,1) = vbtsdta(ii,1) + vsdta(ii,jk,1)*fse3v(ji,jj,jk)
870                      END DO
871                      sshsdta(ii,1) = sshn(ji,jj) * tmask(ii,jj,1)
872                   END DO
873                END DO
874             ENDIF
875          ENDIF
876
877       ENDIF        !       END IF kt == nit000
878
879      !!------------------------------------------------------------------------------------
880      ! 2.      Initialize the time we are at. Does this every time the routine is called,
881      !         excepted when nobc_dta = 0
882      !
883      IF( nobc_dta == 0) THEN
884         itimo = 1
885         zxy   = 0
886      ELSE
887         IF(ntobc == 1) THEN
888            itimo = 1
889         ELSE IF (ntobc == 12) THEN      !   BC are monthly
890            ! we assume we have climatology in that case
891            iman  = 12
892            i15   = nday / 16
893            imois = nmonth + i15 - 1
894            IF( imois == 0 )   imois = iman
895            itimo = imois
896         ELSE
897            IF(lwp) WRITE(numout,*) 'data other than constant or monthly',kt
898            iman  = ntobc
899            itimo = FLOOR( kt*rdt / tcobc(1))
900            isrel=kt*rdt
901         ENDIF
902      ENDIF
903
904      ! 2. Read two records in the file if necessary
905      ! ---------------------------------------------
906
907      IF( nobc_dta == 1 .AND. nlecto == 1 ) THEN
908
909         IF( lp_obc_east ) THEN
910            ! ... Read datafile and set sea surface height and barotropic velocity
911            ! ... initialise the sshedta, ubtedta arrays
912            sshedta(:,0) = sshedta(:,1)
913            ubtedta(:,0) = ubtedta(:,1)
914            CALL iom_open ( 'obceast_TS.nc', id_e )
915            CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,1), ktime=ntobc1 )
916            CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,2), ktime=ntobc2 )
917            IF( lk_dynspg_ts ) THEN
918               CALL iom_get (id_e, jpdom_unknown, 'vossurfh', sshedta(:,3), ktime=ntobc2+1 )
919            ENDIF
920            CALL iom_close ( id_e )
921            !
922            CALL iom_open ( 'obceast_U.nc', id_e )
923            CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,1), ktime=ntobc1 )
924            CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,2), ktime=ntobc2 )
925            IF( lk_dynspg_ts ) THEN
926               CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,3), ktime=ntobc2+1 )
927            ENDIF
928            CALL iom_close ( id_e )
929            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
930            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
931               WRITE(numout,*)
932               WRITE(numout,*) ' Read East OBC barotropic data records ', ntobc1, ntobc2
933               iprint = (jpjef-jpjed+1)/20 +1
934               WRITE(numout,*)
935               WRITE(numout,*) ' Sea surface height record 1'
936               CALL prihre( sshedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout )
937               WRITE(numout,*)
938               WRITE(numout,*) ' Normal transport (m2/s) record 1',iprint
939               CALL prihre( ubtedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout )
940            ENDIF
941         ENDIF
942
943         IF( lp_obc_west ) THEN
944            ! ... Read datafile and set temperature, salinity and normal velocity
945            ! ... initialise the swdta, twdta, uwdta arrays
946            sshwdta(:,0) = sshwdta(:,1)
947            ubtwdta(:,0) = ubtwdta(:,1)
948            CALL iom_open ( 'obcwest_TS.nc', id_w )
949            CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,1), ktime=ntobc1 )
950            CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,2), ktime=ntobc2 )
951            IF( lk_dynspg_ts ) THEN
952               CALL  ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,3), ktime=ntobc2+1 )
953            ENDIF
954            CALL iom_close ( id_w )
955            !
956            CALL iom_open ( 'obcwest_U.nc', id_w )
957            CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,1), ktime=ntobc1 )
958            CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,2), ktime=ntobc2 )
959            IF( lk_dynspg_ts ) THEN
960               CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,3), ktime=ntobc2+1 )
961            ENDIF
962            CALL iom_close ( id_w )
963            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
964            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
965               WRITE(numout,*)
966               WRITE(numout,*) ' Read West OBC barotropic data records ', ntobc1, ntobc2
967               iprint = (jpjwf-jpjwd+1)/20 +1
968               WRITE(numout,*)
969               WRITE(numout,*) ' Sea surface height record 1  - printout surface level'
970               CALL prihre( sshwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout )
971               WRITE(numout,*)
972               WRITE(numout,*) ' Normal transport (m2/s) record 1'
973               CALL prihre( ubtwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout )
974            ENDIF
975         ENDIF
976
977         IF( lp_obc_north) THEN
978            ! ... Read datafile and set sea surface height and barotropic velocity
979            ! ... initialise the sshndta, ubtndta arrays
980            sshndta(:,0) = sshndta(:,1)
981            vbtndta(:,0) = vbtndta(:,1)
982            CALL iom_open ( 'obcnorth_TS.nc', id_n )
983            CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,1), ktime=ntobc1 )
984            CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,2), ktime=ntobc2 )
985            IF( lk_dynspg_ts ) THEN
986               CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,3), ktime=ntobc2+1 )
987            ENDIF
988            CALL iom_close ( id_n )
989
990            CALL iom_open ( 'obcnorth_V.nc', id_n )
991            CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,1), ktime=ntobc1 )
992            CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,2), ktime=ntobc2 )
993            IF( lk_dynspg_ts ) THEN
994               CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,3), ktime=ntobc2+1 )
995            ENDIF
996            CALL iom_close ( id_n )
997
998            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
999            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
1000               WRITE(numout,*)
1001               WRITE(numout,*) ' Read North OBC barotropic data records ', ntobc1, ntobc2
1002               iprint = (jpinf-jpind+1)/20 +1
1003               WRITE(numout,*)
1004               WRITE(numout,*) ' Sea surface height record 1  - printout surface level'
1005               CALL prihre( sshndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout )
1006               WRITE(numout,*)
1007               WRITE(numout,*) ' Normal transport (m2/s) record 1'
1008               CALL prihre( vbtndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout )
1009            ENDIF
1010         ENDIF
1011
1012         IF( lp_obc_south) THEN
1013            ! ... Read datafile and set sea surface height and barotropic velocity
1014            ! ... initialise the sshsdta, ubtsdta arrays
1015            sshsdta(:,0) = sshsdta(:,1)
1016            vbtsdta(:,0) = vbtsdta(:,1)
1017            CALL iom_open ( 'obcsouth_TS.nc', id_s )
1018            CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,1), ktime=ntobc1 )
1019            CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,2), ktime=ntobc2 )
1020            IF( lk_dynspg_ts ) THEN
1021               CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,3), ktime=ntobc2+1 )
1022            ENDIF
1023            CALL iom_close ( id_s )
1024
1025            CALL iom_open ( 'obcsouth_V.nc', id_s )
1026            CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,1), ktime=ntobc1 )
1027            CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,2), ktime=ntobc2 )
1028            IF( lk_dynspg_ts ) THEN
1029               CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,3), ktime=ntobc2+1 )
1030            ENDIF
1031            CALL iom_close ( id_s )
1032               
1033            ! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
1034            IF( lwp .AND.  ( kt == nit000 .OR. nprint /= 0 )  ) THEN
1035               WRITE(numout,*)
1036               WRITE(numout,*) ' Read South OBC barotropic data records ', ntobc1, ntobc2
1037               iprint = (jpisf-jpisd+1)/20 +1
1038               WRITE(numout,*)
1039               WRITE(numout,*) ' Sea surface height record 1  - printout surface level'
1040               CALL prihre( sshsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout )
1041               WRITE(numout,*)
1042               WRITE(numout,*) ' Normal transport (m2/s) record 1'
1043               CALL prihre( vbtsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout )
1044            ENDIF
1045         ENDIF
1046
1047       ENDIF        !      end of the test on the condition to read or not the files
1048
1049      ! 3.  Call at every time step : Linear interpolation of BCs to current time step
1050      ! ----------------------------------------------------------------------
1051
1052       IF( lk_dynspg_ts ) THEN
1053          isrel = (kt-1)*rdt + kbt*rdtbt
1054
1055          IF( nobc_dta == 1 ) THEN
1056             isrel = (kt-1)*rdt + kbt*rdtbt
1057             itimo  = FLOOR(  kt*rdt    / (tcobc(2)-tcobc(1)) )
1058             itimom = FLOOR( (kt-1)*rdt / (tcobc(2)-tcobc(1)) )
1059             itimop = FLOOR( (kt+1)*rdt / (tcobc(2)-tcobc(1)) )
1060             IF( itimom == itimo .AND. itimop == itimo ) THEN
1061                itobcm = ntobc1
1062                itobcp = ntobc2
1063
1064             ELSEIF ( itimom <= itimo .AND. itimop == itimo ) THEN
1065                IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN
1066                   itobcm = ntobc1-1
1067                   itobcp = ntobc2-1
1068                ELSE
1069                   itobcm = ntobc1
1070                   itobcp = ntobc2
1071                ENDIF
1072
1073             ELSEIF ( itimom == itimo .AND. itimop >= itimo ) THEN
1074                IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN
1075                   itobcm = ntobc1
1076                   itobcp = ntobc2
1077                ELSE
1078                   itobcm = ntobc1+1
1079                   itobcp = ntobc2+1
1080                ENDIF
1081
1082             ELSEIF ( itimom == itimo-1 .AND. itimop == itimo+1 ) THEN
1083                IF(  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN
1084                   itobcm = ntobc1-1
1085                   itobcp = ntobc2-1
1086                ELSEIF (  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN
1087                   itobcm = ntobc1
1088                   itobcp = ntobc2
1089                ELSEIF (  FLOOR( isrel / (tcobc(2)-tcobc(1)) ) == itimop ) THEN
1090                   itobcm = ntobc1+1
1091                   itobcp = ntobc2+2
1092                ELSE
1093                   IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 1?'
1094                ENDIF
1095             ELSE
1096                IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 2?'
1097             ENDIF
1098
1099          ENDIF
1100
1101       ELSE IF( lk_dynspg_exp ) THEN
1102          isrel=kt*rdt
1103          itobcm = ntobc1
1104          itobcp = ntobc2
1105       ENDIF
1106
1107       IF( ntobc == 1 .OR. nobc_dta == 0 ) THEN
1108          zxy = 0.e0
1109       ELSE IF( ntobc == 12 ) THEN
1110          zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
1111       ELSE
1112          zxy = (tcobc(itobcm)-FLOAT(isrel)) / (tcobc(itobcm)-tcobc(itobcp))
1113       ENDIF
1114
1115       IF( lp_obc_east ) THEN           !  fills sshfoe, ubtfoe (local to each processor)
1116          DO jj = nje0p1, nje1m1
1117             ij = jj -1 + njmpp
1118             sshfoe(jj) = ( zxy * sshedta(ij,2) + (1.-zxy) * sshedta(ij,1) ) * temsk(jj,1)
1119             ubtfoe(jj) = ( zxy * ubtedta(ij,2) + (1.-zxy) * ubtedta(ij,1) ) * uemsk(jj,1)
1120          END DO
1121       ENDIF
1122
1123       IF( lp_obc_west) THEN            !  fills sshfow, ubtfow (local to each processor)
1124          DO jj = njw0p1, njw1m1
1125             ij = jj -1 + njmpp
1126             sshfow(jj) = ( zxy * sshwdta(ij,2) + (1.-zxy) * sshwdta(ij,1) ) * twmsk(jj,1)
1127             ubtfow(jj) = ( zxy * ubtwdta(ij,2) + (1.-zxy) * ubtwdta(ij,1) ) * uwmsk(jj,1)
1128          END DO
1129       ENDIF
1130
1131       IF( lp_obc_north) THEN           !  fills sshfon, vbtfon (local to each processor)
1132          DO ji = nin0p1, nin1m1
1133             ii = ji -1 + nimpp
1134             sshfon(ji) = ( zxy * sshndta(ii,2) + (1.-zxy) * sshndta(ii,1) ) * tnmsk(ji,1)
1135             vbtfon(ji) = ( zxy * vbtndta(ii,2) + (1.-zxy) * vbtndta(ii,1) ) * vnmsk(ji,1)
1136          END DO
1137       ENDIF
1138
1139       IF( lp_obc_south) THEN           !  fills sshfos, vbtfos (local to each processor)
1140          DO ji = nis0p1, nis1m1
1141             ii = ji -1 + nimpp
1142             sshfos(ji) = ( zxy * sshsdta(ii,2) + (1.-zxy) * sshsdta(ii,1) ) * tsmsk(ji,1)
1143             vbtfos(ji) = ( zxy * vbtsdta(ii,2) + (1.-zxy) * vbtsdta(ii,1) ) * vsmsk(ji,1)
1144          END DO
1145       ENDIF
1146
1147   END SUBROUTINE obc_dta_bt
1148
1149#else
1150   !!-----------------------------------------------------------------------------
1151   !!   Default option
1152   !!-----------------------------------------------------------------------------
1153   SUBROUTINE obc_dta_bt ( kt, kbt )       ! Empty routine
1154      !! * Arguments
1155      INTEGER,INTENT(in) :: kt
1156      INTEGER, INTENT( in ) ::   kbt         ! barotropic ocean time-step index
1157      WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt
1158      WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kbt
1159   END SUBROUTINE obc_dta_bt
1160#endif
1161
1162#else
1163   !!--------------------------------------------------------------------
1164   !!  default option  :  Dummy module    NO Open Boundary Conditions
1165   !!--------------------------------------------------------------------
1166CONTAINS
1167   SUBROUTINE obc_dta( kt )             ! Dummy routine
1168     INTEGER, INTENT (in) :: kt
1169     WRITE(*,*) 'obc_dta: You should not have seen this print! error?', kt
1170   END SUBROUTINE obc_dta
1171   SUBROUTINE obc_dta_bt( kt, jn)             ! Dummy routine
1172     INTEGER, INTENT (in) :: kt, jn
1173     WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', kt
1174     WRITE(*,*) 'obc_dta_bt: You should not have seen this print! error?', jn
1175   END SUBROUTINE obc_dta_bt
1176#endif
1177
1178   !!=====================================================================
1179END MODULE obcdta
Note: See TracBrowser for help on using the repository browser.