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

Context Navigation

source: trunk/NEMO/OPA_SRC/OBC/obcdta.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: 52.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: