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 @ 1125

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

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

Download in other formats: