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

Last change on this file since 1151 was 1151, checked in by rblod, 16 years ago
Add new OBC, see ticket #224
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 64.2 KB

Line
1	MODULE obcdta
2	!!==============================================================================
3	!! * MODULE obcdta *
4	!! Open boundary data : read the data for the open boundaries.
5	!!==============================================================================
6	#if defined key_obc
7	!!------------------------------------------------------------------------------
8	!! 'key_obc' : Open Boundary Conditions
9	!!------------------------------------------------------------------------------
10	!! obc_dta : read u, v, t, s data along each open boundary
11	!! obc_dta_psi : read psi data along each open boundary (rigid lid only)
12	!!------------------------------------------------------------------------------
13	!! * Modules used
14	USE oce ! ocean dynamics and tracers
15	USE dom_oce ! ocean space and time domain
16	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
17	USE phycst ! physical constants
18	USE obc_oce ! ocean open boundary conditions
19	USE daymod ! calendar
20	USE in_out_manager ! I/O logical units
21	USE lib_mpp ! distributed memory computing
22	USE dynspg_oce
23	USE ioipsl ! now only for ymds2ju function
24	USE iom !
25	USE obccli ! used when dynspg_rl
26
27	IMPLICIT NONE
28	PRIVATE
29
30	!! * Accessibility
31	PUBLIC obc_dta ! routines called by step.F90
32	PUBLIC obc_dta_bt ! routines called by dynspg_ts.F90
33
34	!! * Shared module variables
35	REAL(wp), DIMENSION(2) :: zjcnes_obc !
36	REAL(wp), DIMENSION(:), ALLOCATABLE :: ztcobc
37	REAL(wp) :: rdt_obc
38	REAL(wp) :: zjcnes
39	INTEGER :: imm0, iyy0, idd0, iyy, imm, idd
40	INTEGER :: nt_a=2, nt_b=1, itobc, ndate0_cnes, nday_year0
41	INTEGER :: itobce, itobcw, itobcs, itobcn, itobc_b ! number of time steps in OBC files
42	REAL(wp) :: zraaj=365.
43
44	INTEGER :: &
45	ntobc , & !: where we are in the obc file
46	ntobc_b , & !: first record used
47	ntobc_a !: second record used
48
49	CHARACTER (len=40) :: & ! name of data files
50	cl_obc_eTS , cl_obc_eU, &
51	cl_obc_wTS , cl_obc_wU, &
52	cl_obc_nTS , cl_obc_nV, &
53	cl_obc_sTS , cl_obc_sV
54
55	! arrays used for interpolating time dependent data on the boundaries
56	REAL(wp), DIMENSION(jpj,jpk,jptobc) :: uedta, vedta, tedta, sedta ! East
57	REAL(wp), DIMENSION(jpj,jpk,jptobc) :: uwdta, vwdta, twdta, swdta ! West
58	REAL(wp), DIMENSION(jpi,jpk,jptobc) :: undta, vndta, tndta, sndta ! North
59	REAL(wp), DIMENSION(jpi,jpk,jptobc) :: usdta, vsdta, tsdta, ssdta ! South
60
61	LOGICAL, DIMENSION (jpj,jpk ) :: ltemsk=.TRUE., luemsk=.TRUE., lvemsk=.TRUE. ! boolean msks
62	LOGICAL, DIMENSION (jpj,jpk ) :: ltwmsk=.TRUE., luwmsk=.TRUE., lvwmsk=.TRUE. ! used for outliers
63	LOGICAL, DIMENSION (jpi,jpk ) :: ltnmsk=.TRUE., lunmsk=.TRUE., lvnmsk=.TRUE. ! checks
64	LOGICAL, DIMENSION (jpi,jpk ) :: ltsmsk=.TRUE., lusmsk=.TRUE., lvsmsk=.TRUE.
65
66	!! * Substitutions
67	# include "obc_vectopt_loop_substitute.h90"
68	!!---------------------------------------------------------------------------------
69	!! OPA 9.0 , LODYC-IPSL (2003)
70	!!---------------------------------------------------------------------------------
71
72	CONTAINS
73
74	SUBROUTINE obc_dta( kt )
75	!!---------------------------------------------------------------------------
76	!! * SUBROUTINE obc_dta *
77	!!
78	!! ** Purpose : Find the climatological boundary arrays for the specified date,
79	!! The boundary arrays are netcdf files. Three possible cases:
80	!! - one time frame only in the file (time dimension = 1).
81	!! in that case the boundary data does not change in time.
82	!! - many time frames. In that case, if we have 12 frames
83	!! we assume monthly fields.
84	!! Else, we assume that time_counter is in seconds
85	!! since the beginning of either the current year or a reference
86	!! year given in the namelist.
87	!! (no check is done so far but one would have to check the "unit"
88	!! attribute of variable time_counter).
89	!!
90	!!
91	!! History :
92	!! ! 98-05 (J.M. Molines) Original code
93	!! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
94	!!
95	!! 9.0 ! 04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input
96	!! ! 2007-2008 (C. Langlais, P. Mathiot, J.M. Molines) high frequency boundaries data
97	!!---------------------------------------------------------------------------
98	!! * Arguments
99	INTEGER, INTENT( in ) :: kt ! ocean time-step index
100
101	!! * Local declarations
102	INTEGER :: ji, jj, jk, ii, ij ,it ! dummy loop indices
103	INTEGER :: ikprint ! frequency for printouts.
104	INTEGER, SAVE :: immfile, iyyfile !
105	INTEGER :: nt ! record indices (incrementation)
106	INTEGER :: istop ! local error check
107
108	REAL(wp) :: zxy, znum, zden ! time interpolation weight
109
110	! variables for the julian day calculation
111	INTEGER :: iyear, imonth, iday
112	REAL(wp) :: zsec , zjulian, zjuliancnes
113
114	! IOM STUFF
115	INTEGER :: idvar, id_e, id_w, id_n, id_s, id_x ! file identifiers
116	INTEGER, DIMENSION(1) :: itmp
117	CHARACTER(LEN=25) :: cl_vname
118
119	!!---------------------------------------------------------------------------
120	IF( lk_dynspg_rl ) THEN
121	CALL obc_dta_psi( kt ) ! update bsf data at open boundaries
122	IF ( nobc_dta == 1 .AND. kt == nit000 ) CALL ctl_stop( 'obcdta : time-variable psi boundary data not allowed yet' )
123	ENDIF
124
125	! 0. initialisation :
126	! --------------------
127	IF ( kt == nit000 ) CALL obc_dta_ini ( kt )
128	IF ( nobc_dta == 0 ) RETURN ! already done in obc_dta_ini
129	IF ( itobc == 1 ) RETURN ! case of only one time frame in file done in obc_dta_ini
130
131	! in the following code, we assume that obc data are read from files, with more than 1 time frame in it
132
133	iyyfile=iyy ; immfile = 00 ! set component of the current file name
134	IF ( cffile /= 'annual') immfile = imm !
135	IF ( ln_obc_clim ) iyyfile = 0000 ! assume that climatological files are labeled y0000
136
137	! 1. Synchronize time of run with time of data files
138	!---------------------------------------------------
139	! nday_year is the day number in the current year ( 1 for 01/01 )
140	zsec=MOD( (kt-nit000)rdt - (nday_year - nday_year0 )rday, rday ) ! number of seconds in the current day
141	IF (ln_obc_clim) THEN
142	zjcnes = nday_year - 1 + zsec/rday
143	ELSE
144	zjcnes = zjcnes + rdt/rday
145	ENDIF
146
147	! look for 'before' record number in the current file
148	ntobc = nrecbef () ! this function return the record number for 'before', relative to zjcnes
149
150	IF (MOD(kt-1,10)==0) THEN
151	IF (lwp) WRITE(numout,*) 'kt= ',kt,' zjcnes =', zjcnes,' ndastp =',ndastp, 'mm =',imm
152	END IF
153
154	! 2. read a new data if necessary
155	!--------------------------------
156	IF ( ntobc /= ntobc_b ) THEN
157	! we need to read the 'after' record
158	! swap working index:
159	nt=nt_b ; nt_b=nt_a ; nt_a=nt
160	ntobc_b = ntobc
161
162	! new record number :
163	ntobc_a = ntobc_a + 1
164
165	! all tricky things related to record number, changing files etc... are managed by obc_read
166
167	CALL obc_read (kt, nt_a, ntobc_a, iyyfile, immfile )
168
169	! update zjcnes_obc
170	zjcnes_obc(nt_b)= ztcobc(ntobc_b)
171	zjcnes_obc(nt_a)= ztcobc(ntobc_a)
172	ENDIF
173
174	! 3. interpolation at each time step
175	! ------------------------------------
176	IF ( ln_obc_clim) THEN
177	znum= MOD(zjcnes - zjcnes_obc(nt_b), zraaj ) ; IF ( znum < 0 ) znum = znum + zraaj
178	zden= MOD(zjcnes_obc(nt_a) - zjcnes_obc(nt_b), zraaj ) ; IF ( zden < 0 ) zden = zden + zraaj
179	ELSE
180	znum= zjcnes - zjcnes_obc(nt_b)
181	zden= zjcnes_obc(nt_a) - zjcnes_obc(nt_b)
182	ENDIF
183	zxy = znum / zden
184
185	IF( lp_obc_east ) THEN
186	! fills sfoe, tfoe, ufoe ,vfoe
187	sfoe(:,:) = zxy * sedta (:,:,nt_a) + (1. - zxy)*sedta(:,:,nt_b)
188	tfoe(:,:) = zxy * tedta (:,:,nt_a) + (1. - zxy)*tedta(:,:,nt_b)
189	ufoe(:,:) = zxy * uedta (:,:,nt_a) + (1. - zxy)*uedta(:,:,nt_b)
190	vfoe(:,:) = zxy * vedta (:,:,nt_a) + (1. - zxy)*vedta(:,:,nt_b)
191	ENDIF
192
193	IF( lp_obc_west) THEN
194	! fills sfow, tfow, ufow ,vfow
195	sfow(:,:) = zxy * swdta (:,:,nt_a) + (1. - zxy)*swdta(:,:,nt_b)
196	tfow(:,:) = zxy * twdta (:,:,nt_a) + (1. - zxy)*twdta(:,:,nt_b)
197	ufow(:,:) = zxy * uwdta (:,:,nt_a) + (1. - zxy)*uwdta(:,:,nt_b)
198	vfow(:,:) = zxy * vwdta (:,:,nt_a) + (1. - zxy)*vwdta(:,:,nt_b)
199	ENDIF
200
201	IF( lp_obc_north) THEN
202	! fills sfon, tfon, ufon ,vfon
203	sfon(:,:) = zxy * sndta (:,:,nt_a) + (1. - zxy)*sndta(:,:,nt_b)
204	tfon(:,:) = zxy * tndta (:,:,nt_a) + (1. - zxy)*tndta(:,:,nt_b)
205	ufon(:,:) = zxy * undta (:,:,nt_a) + (1. - zxy)*undta(:,:,nt_b)
206	vfon(:,:) = zxy * vndta (:,:,nt_a) + (1. - zxy)*vndta(:,:,nt_b)
207	ENDIF
208
209	IF( lp_obc_south) THEN
210	! fills sfos, tfos, ufos ,vfos
211	sfos(:,:) = zxy * ssdta (:,:,nt_a) + (1. - zxy)*ssdta(:,:,nt_b)
212	tfos(:,:) = zxy * tsdta (:,:,nt_a) + (1. - zxy)*tsdta(:,:,nt_b)
213	ufos(:,:) = zxy * usdta (:,:,nt_a) + (1. - zxy)*usdta(:,:,nt_b)
214	vfos(:,:) = zxy * vsdta (:,:,nt_a) + (1. - zxy)*vsdta(:,:,nt_b)
215	ENDIF
216	END SUBROUTINE obc_dta
217
218
219	SUBROUTINE obc_dta_ini (kt)
220	!!-----------------------------------------------------------------------------
221	!! * SUBROUTINE obc_dta_ini *
222	!!
223	!! ** Purpose :
224	!! When obc_dta first call, realize some data initialization
225	!!
226	!! ** Method :
227	!!
228	!! History :
229	!! 9.0 ! 07-10 (J.M. Molines )
230	!!----------------------------------------------------------------------------
231	!! * Argument
232	INTEGER, INTENT(in) :: kt ! ocean time-step index
233
234	!! * Local declarations
235	INTEGER :: ji,jj, it ! dummy loop indices
236
237	REAL(wp) :: zxy ! time interpolation weight
238
239	INTEGER :: ikprint ! frequency for printouts.
240
241	INTEGER, SAVE :: immfile, iyyfile !
242	INTEGER :: nt ! record indices (incrementation)
243	INTEGER :: istop ! local error check
244
245	! variables for the julian day calculation
246	INTEGER :: iyear, imonth, iday
247	REAL(wp) :: zsec , zjulian, zjuliancnes
248
249	! IOM STUFF
250	INTEGER :: idvar, id_e, id_w, id_n, id_s, id_x ! file identifiers
251	INTEGER, DIMENSION(1) :: itmp
252	CHARACTER(LEN=25) :: cl_vname
253
254	IF(lwp) WRITE(numout,*)
255	IF(lwp) WRITE(numout,*) 'obc_dta : find boundary data'
256	IF(lwp) WRITE(numout,*) '~~~~~~~'
257	IF (lwp) THEN
258	IF ( nobc_dta == 0 ) THEN
259	WRITE(numout,*) ' OBC data taken from initial conditions.'
260	ELSE
261	WRITE(numout,*) ' OBC data taken from netcdf files.'
262	ENDIF
263	ENDIF
264	nday_year0 = nday_year ! to remember the day when kt=nit000
265	SELECT CASE (nleapy)
266	CASE ( 0 ) ; zraaj = 365.e0
267	CASE ( 1 ) ; zraaj = 366.e0 !! ERROR TO CORRECT: NOT EVERY YEAR IS LEAP YEAR !! JMM
268	CASE DEFAULT; zraaj = 12.* nleapy
269	END SELECT
270
271	sedta(:,:,:) = 0.e0 ; tedta(:,:,:) = 0.e0 ; uedta(:,:,:) = 0.e0 ; vedta(:,:,:) = 0.e0 ! East
272	swdta(:,:,:) = 0.e0 ; twdta(:,:,:) = 0.e0 ; uwdta(:,:,:) = 0.e0 ; vwdta(:,:,:) = 0.e0 ! West
273	sndta(:,:,:) = 0.e0 ; tndta(:,:,:) = 0.e0 ; undta(:,:,:) = 0.e0 ; vndta(:,:,:) = 0.e0 ! North
274	ssdta(:,:,:) = 0.e0 ; tsdta(:,:,:) = 0.e0 ; usdta(:,:,:) = 0.e0 ; vsdta(:,:,:) = 0.e0 ! South
275
276	sfoe(:,:) = 0.e0 ; tfoe(:,:) = 0.e0 ; ufoe(:,:) = 0.e0 ; vfoe(:,:) = 0.e0 ! East
277	sfow(:,:) = 0.e0 ; tfow(:,:) = 0.e0 ; ufow(:,:) = 0.e0 ; vfow(:,:) = 0.e0 ! West
278	sfon(:,:) = 0.e0 ; tfon(:,:) = 0.e0 ; ufon(:,:) = 0.e0 ; vfon(:,:) = 0.e0 ! North
279	sfos(:,:) = 0.e0 ; tfos(:,:) = 0.e0 ; ufos(:,:) = 0.e0 ; vfos(:,:) = 0.e0 ! South
280
281	IF (nobc_dta == 0 ) THEN ! boundary data are the initial data of this run (set only at nit000)
282	IF (lp_obc_east) THEN ! East
283	DO ji = nie0 , nie1
284	sfoe(nje0p1:nje1m1,:) = temsk(nje0p1:nje1m1,:) * sn (ji+1 , nje0p1:nje1m1 , :)
285	tfoe(nje0p1:nje1m1,:) = temsk(nje0p1:nje1m1,:) * tn (ji+1 , nje0p1:nje1m1 , :)
286	ufoe(nje0p1:nje1m1,:) = uemsk(nje0p1:nje1m1,:) * un (ji , nje0p1:nje1m1 , :)
287	vfoe(nje0p1:nje1m1,:) = vemsk(nje0p1:nje1m1,:) * vn (ji+1 , nje0p1:nje1m1 , :)
288	END DO
289	ENDIF
290
291	IF (lp_obc_west) THEN ! West
292	DO ji = niw0 , niw1
293	sfow(njw0p1:njw1m1,:) = twmsk(njw0p1:njw1m1,:) * sn (ji , njw0p1:njw1m1 , :)
294	tfow(njw0p1:njw1m1,:) = twmsk(njw0p1:njw1m1,:) * tn (ji , njw0p1:njw1m1 , :)
295	ufow(njw0p1:njw1m1,:) = uwmsk(njw0p1:njw1m1,:) * un (ji , njw0p1:njw1m1 , :)
296	vfow(njw0p1:njw1m1,:) = vwmsk(njw0p1:njw1m1,:) * vn (ji , njw0p1:njw1m1 , :)
297	END DO
298	ENDIF
299
300	IF (lp_obc_north) THEN ! North
301	DO jj = njn0 , njn1
302	sfon(nin0p1:nin1m1,:) = tnmsk(nin0p1:nin1m1,:) * sn (nin0p1:nin1m1 , jj+1 , :)
303	tfon(nin0p1:nin1m1,:) = tnmsk(nin0p1:nin1m1,:) * tn (nin0p1:nin1m1 , jj+1 , :)
304	ufon(nin0p1:nin1m1,:) = unmsk(nin0p1:nin1m1,:) * un (nin0p1:nin1m1 , jj+1 , :)
305	vfon(nin0p1:nin1m1,:) = vnmsk(nin0p1:nin1m1,:) * vn (nin0p1:nin1m1 , jj , :)
306	END DO
307	ENDIF
308
309	IF (lp_obc_south) THEN ! South
310	DO jj = njs0 , njs1
311	sfos(nis0p1:nis1m1,:) = tsmsk(nis0p1:nis1m1,:) * sn (nis0p1:nis1m1 , jj , :)
312	tfos(nis0p1:nis1m1,:) = tsmsk(nis0p1:nis1m1,:) * tn (nis0p1:nis1m1 , jj , :)
313	ufos(nis0p1:nis1m1,:) = usmsk(nis0p1:nis1m1,:) * un (nis0p1:nis1m1 , jj , :)
314	vfos(nis0p1:nis1m1,:) = vsmsk(nis0p1:nis1m1,:) * vn (nis0p1:nis1m1 , jj , :)
315	END DO
316	ENDIF
317	RETURN ! exit the routine all is done
318	ENDIF ! nobc_dta = 0
319
320	!!!! In the following OBC data are read from files.
321	! all logical-mask are initialzed to true when declared
322	WHERE ( temsk == 0 ) ltemsk=.FALSE.
323	WHERE ( uemsk == 0 ) luemsk=.FALSE.
324	WHERE ( vemsk == 0 ) lvemsk=.FALSE.
325
326	WHERE ( twmsk == 0 ) ltwmsk=.FALSE.
327	WHERE ( uwmsk == 0 ) luwmsk=.FALSE.
328	WHERE ( vwmsk == 0 ) lvwmsk=.FALSE.
329
330	WHERE ( tnmsk == 0 ) ltnmsk=.FALSE.
331	WHERE ( unmsk == 0 ) lunmsk=.FALSE.
332	WHERE ( vnmsk == 0 ) lvnmsk=.FALSE.
333
334	WHERE ( tsmsk == 0 ) ltsmsk=.FALSE.
335	WHERE ( usmsk == 0 ) lusmsk=.FALSE.
336	WHERE ( vsmsk == 0 ) lvsmsk=.FALSE.
337
338	iyear=1950; imonth=01; iday=01; zsec=0.
339	! zjuliancnes : julian day corresonding to 01/01/1950
340	CALL ymds2ju(iyear, imonth, iday,zsec , zjuliancnes)
341
342	!current year and curent month
343	iyy=INT(ndastp/10000) ; imm=INT((ndastp -iyy10000)/100) ; idd=(ndastp-iyy10000-imm*100)
344	IF (iyy < 1900) iyy = iyy+1900 ! always assume that years are on 4 digits.
345	CALL ymds2ju(iyy, imm, idd ,zsec , zjulian)
346	ndate0_cnes = zjulian - zjuliancnes ! jcnes day when call to obc_dta_ini
347
348	iyyfile=iyy ; immfile=0 ! set component of the current file name
349	IF ( cffile /= 'annual') immfile=imm
350	IF ( ln_obc_clim) iyyfile = 0 ! assume that climatological files are labeled y0000
351
352	CALL obc_dta_chktime ( iyyfile, immfile )
353
354	IF ( itobc == 1 ) THEN
355	! in this case we will provide boundary data only once.
356	nt_a=1 ; ntobc_a=1
357	CALL obc_read (nit000, nt_a, ntobc_a, iyyfile, immfile)
358	IF( lp_obc_east ) THEN
359	! fills sfoe, tfoe, ufoe ,vfoe
360	sfoe(:,:) = sedta (:,:,1) ; tfoe(:,:) = tedta (:,:,1)
361	ufoe(:,:) = uedta (:,:,1) ; vfoe(:,:) = vedta (:,:,1)
362	ENDIF
363
364	IF( lp_obc_west) THEN
365	! fills sfow, tfow, ufow ,vfow
366	sfow(:,:) = swdta (:,:,1) ; tfow(:,:) = twdta (:,:,1)
367	ufow(:,:) = uwdta (:,:,1) ; vfow(:,:) = vwdta (:,:,1)
368	ENDIF
369
370	IF( lp_obc_north) THEN
371	! fills sfon, tfon, ufon ,vfon
372	sfon(:,:) = sndta (:,:,1) ; tfon(:,:) = tndta (:,:,1)
373	ufon(:,:) = undta (:,:,1) ; vfon(:,:) = vndta (:,:,1)
374	ENDIF
375
376	IF( lp_obc_south) THEN
377	! fills sfos, tfos, ufos ,vfos
378	sfos(:,:) = ssdta (:,:,1) ; tfos(:,:) = tsdta (:,:,1)
379	ufos(:,:) = usdta (:,:,1) ; vfos(:,:) = vsdta (:,:,1)
380	ENDIF
381	RETURN ! we go out of obc_dta_ini -------------------------------------->>>>>
382	ENDIF
383
384	! nday_year is the day number in the current year ( 1 for 01/01 )
385	! we suppose that we always start from the begining of a day
386	! zsec=MOD( (kt-nit000)rdt - (nday_year - nday_year0 )rday, rday ) ! number of seconds in the current day
387	zsec=0.e0 ! here, kt=nit000, nday_year = ndat_year0
388
389	IF (ln_obc_clim) THEN
390	zjcnes = nday_year - 1 + zsec/rday ! for clim file time is in days in a year
391	ELSE
392	zjcnes = ndate0_cnes + (nday_year - nday_year0 ) + zsec/rday
393	ENDIF
394
395	! look for 'before' record number in the current file
396	ntobc = nrecbef ()
397
398	IF (lwp) WRITE(numout,*) 'obc files frequency :',cffile
399	IF (lwp) WRITE(numout,*) ' zjcnes0 =',zjcnes,' ndastp0 =',ndastp
400	IF (lwp) WRITE(numout,*) ' annee0 ',iyy,' month0 ', imm,' day0 ', idd
401	IF (lwp) WRITE(numout,*) 'first file open :',cl_obc_nTS
402
403	! record initialisation
404	!--------------------
405	nt_b = 1 ; nt_a = 2
406
407	ntobc_a = ntobc + 1
408	ntobc_b = ntobc
409
410	CALL obc_read (kt, nt_b, ntobc_b, iyyfile, immfile) ! read 'before' fields
411	CALL obc_read (kt, nt_a, ntobc_a, iyyfile, immfile) ! read 'after' fields
412
413	zjcnes_obc(nt_b)= ztcobc(ntobc_b)
414	zjcnes_obc(nt_a)= ztcobc(ntobc_a)
415	!
416	END SUBROUTINE obc_dta_ini
417
418
419	SUBROUTINE obc_dta_chktime (kyyfile, kmmfile)
420	!
421	! check the number of time steps in the files and read ztcobc
422	!
423	! * Arguments
424	INTEGER, INTENT(in) :: kyyfile, kmmfile
425	! * local variables
426	INTEGER :: istop ! error control
427	INTEGER :: ji ! dummy loop index
428
429	INTEGER :: idvar, id_e, id_w, id_n, id_s, id_x ! file identifiers
430	INTEGER, DIMENSION(1) :: itmp
431	CHARACTER(LEN=25) :: cl_vname
432
433	ntobc_a = 0; itobce =0 ; itobcw = 0; itobcn = 0; itobcs = 0
434	! build file name
435	WRITE(cl_obc_eTS ,'("obc_east_TS_y",i4.4,"m",i2.2,".nc")' ) kyyfile,kmmfile
436	WRITE(cl_obc_wTS ,'("obc_west_TS_y",i4.4,"m",i2.2,".nc")' ) kyyfile,kmmfile
437	WRITE(cl_obc_nTS ,'("obc_north_TS_y",i4.4,"m",i2.2,".nc")' ) kyyfile,kmmfile
438	WRITE(cl_obc_sTS ,'("obc_south_TS_y",i4.4,"m",i2.2,".nc")' ) kyyfile,kmmfile
439
440	cl_vname = 'time_counter'
441	IF ( lp_obc_east ) THEN
442	CALL iom_open ( cl_obc_eTS , id_e )
443	idvar = iom_varid( id_e, cl_vname, kdimsz = itmp ); itobce=itmp(1)
444	ENDIF
445	IF ( lp_obc_west ) THEN
446	CALL iom_open ( cl_obc_wTS , id_w )
447	idvar = iom_varid( id_w, cl_vname, kdimsz = itmp ) ; itobcw=itmp(1)
448	ENDIF
449	IF ( lp_obc_north ) THEN
450	CALL iom_open ( cl_obc_nTS , id_n )
451	idvar = iom_varid( id_n, cl_vname, kdimsz = itmp ) ; itobcn=itmp(1)
452	ENDIF
453	IF ( lp_obc_south ) THEN
454	CALL iom_open ( cl_obc_sTS , id_s )
455	idvar = iom_varid( id_s, cl_vname, kdimsz = itmp ) ; itobcs=itmp(1)
456	ENDIF
457
458	itobc = MAX( itobce, itobcw, itobcn, itobcs )
459	istop = 0
460	IF ( lp_obc_east .AND. itobce /= itobc ) istop = istop+1
461	IF ( lp_obc_west .AND. itobcw /= itobc ) istop = istop+1
462	IF ( lp_obc_north .AND. itobcn /= itobc ) istop = istop+1
463	IF ( lp_obc_south .AND. itobcs /= itobc ) istop = istop+1
464	nstop = nstop + istop
465
466	IF ( istop /= 0 ) THEN
467	WRITE(ctmp1,*) ' east, west, north, south: ', itobce, itobcw, itobcn, itobcs
468	CALL ctl_stop( 'obcdta : all files must have the same number of time steps', ctmp1 )
469	ENDIF
470
471	IF ( itobc == 1 ) THEN
472	IF (lwp) THEN
473	WRITE(numout,*) ' obcdta found one time step only in the OBC files'
474	IF (ln_obc_clim) THEN
475	! OK no problem
476	ELSE
477	ln_obc_clim=.true.
478	WRITE(numout,*) ' we force ln_obc_clim to T'
479	ENDIF
480	ENDIF
481	ELSE
482	IF ( ALLOCATED(ztcobc) ) DEALLOCATE ( ztcobc )
483	ALLOCATE (ztcobc(itobc))
484	DO ji=1,1 ! use a dummy loop to read ztcobc only once
485	IF ( lp_obc_east ) THEN
486	CALL iom_gettime ( id_e, ztcobc, cl_vname ) ; CALL iom_close (id_e) ; EXIT
487	ENDIF
488	IF ( lp_obc_west ) THEN
489	CALL iom_gettime ( id_w, ztcobc, cl_vname ) ; CALL iom_close (id_w) ; EXIT
490	ENDIF
491	IF ( lp_obc_north ) THEN
492	CALL iom_gettime ( id_n, ztcobc, cl_vname ) ; CALL iom_close (id_n) ; EXIT
493	ENDIF
494	IF ( lp_obc_south ) THEN
495	CALL iom_gettime ( id_s, ztcobc, cl_vname ) ; CALL iom_close (id_s) ; EXIT
496	ENDIF
497	END DO
498	rdt_obc = ztcobc(2)-ztcobc(1) ! just an information, not used for any computation
499	IF (lwp) WRITE(numout,*) ' obcdta found', itobc,' time steps in the OBC files'
500	IF (lwp) WRITE(numout,*) ' time step of obc data :', rdt_obc,' days'
501	ENDIF
502	zjcnes = zjcnes - rdt/rday ! trick : zcnes is always incremented by rdt/rday in obc_dta!
503	END SUBROUTINE obc_dta_chktime
504
505	# if defined key_dynspg_rl
506	!!-----------------------------------------------------------------------------
507	!! Rigid-lid
508	!!-----------------------------------------------------------------------------
509
510	SUBROUTINE obc_dta_psi ( kt )
511	!!-----------------------------------------------------------------------------
512	!! * SUBROUTINE obc_dta_psi *
513	!!
514	!! ** Purpose :
515	!! Update the climatological streamfunction OBC at each time step.
516	!! Depends on the user's configuration. Here data are read only once
517	!! at the beginning of the run.
518	!!
519	!! ** Method :
520	!! 1. initialization
521	!! kbsfstart: number of time steps over which increase bsf
522	!! during initialization. This is provided for a smooth start
523	!! in cases where the transport is large (like on the Antarctic
524	!! continent). also note that when kbfstart=1, the transport
525	!! increases a lot in one time step and the precision usually
526	!! required for the solver may not be enough.
527	!! 2. set the time evolution of the climatological barotropic streamfunction
528	!! along the isolated coastlines ( gcbic(jnic) ).
529	!! 3. set the climatological barotropic streamfunction at the boundary.
530	!!
531	!! The last two steps are done only at first step (nit000) or if kt <= kbfstart
532	!!
533	!! History :
534	!! ! 97-08 (G. Madec, J.M. Molines)
535	!! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
536	!!----------------------------------------------------------------------------
537	!! * Arguments
538	INTEGER, INTENT( in ) :: kt ! ocean time-step index
539
540	!! * Local declarations
541	INTEGER :: ji, jj, jnic, jip ! dummy loop indices
542	INTEGER :: inum = 11 ! temporary logical unit
543	INTEGER :: ip, ii, ij, iii, ijj
544	INTEGER :: kbsfstart
545	REAL(wp) :: zsver1, zsver2, zsver3, z2dtr, zcoef
546	!!----------------------------------------------------------------------------
547
548	! 1. initialisation
549	! -----------------
550	kbsfstart = 1
551	zsver1 = bsfic0(1)
552	zsver2 = zsver1
553	IF( kt <= kbsfstart ) THEN
554	zcoef = float(kt)/float(kbsfstart)
555	ELSE
556	zcoef = 1.
557	END IF
558	bsfic(1) = zsver1*zcoef
559	IF( lwp .AND. ( kt <= kbsfstart ) ) THEN
560	IF(lwp) WRITE(numout,*)' '
561	IF(lwp) WRITE(numout,*)'obcdta: spinup phase in obc_dta_psi routine'
562	IF(lwp) WRITE(numout,*)'~~~~~~ it=',kt,' OBC: spinup coef: ', &
563	zcoef, ' and transport: ',bsfic(1)
564	END IF
565
566	zsver2 = bsfic(1)-bsfic(2)
567	zsver3 = bsfic(2)
568
569	! 2. Right hand side of the barotropic elliptic equation (isolated coastlines)
570	! ----------------------------------------------------------------------------
571
572	IF( ( neuler == 0 ) .AND. ( kt == nit000 ) ) THEN
573	z2dtr = 1./rdt
574	ELSE
575	z2dtr = 1./2./rdt
576	END IF
577	! ... bsfb(ii,ij) should be constant but due to the Asselin filter it
578	! ... converges asymptotically towards bsfic(jnic)
579	! ... However, bsfb(ii,ij) is constant along the same coastlines
580	! ... ---> can be improved using an extra array for storing bsficb (before)
581	IF( nbobc > 1 ) THEN
582	DO jnic = 1,nbobc - 1
583	gcbic(jnic) = 0.
584	ip=mnic(0,jnic)
585	DO jip = 1,ip
586	ii = miic(jip,0,jnic)
587	ij = mjic(jip,0,jnic)
588	IF( ii >= nldi+ nimpp - 1 .AND. ii <= nlei+ nimpp - 1 .AND. &
589	ij >= nldj+ njmpp - 1 .AND. ij <= nlej+ njmpp - 1 ) THEN
590	iii=ii-nimpp+1
591	ijj=ij-njmpp+1
592	gcbic(jnic) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
593	END IF
594	END DO
595	END DO
596	END IF
597
598	IF( lk_mpp ) CALL mpp_isl( gcbic, 3 )
599
600	! 3. Update the climatological barotropic function at the boundary
601	! ----------------------------------------------------------------
602
603	IF( lpeastobc ) THEN
604
605	IF( kt == nit000 .OR. kt <= kbsfstart ) THEN
606	OPEN(inum,file='obceastbsf.dta')
607	READ(inum,*)
608	READ(inum,*)
609	READ(inum,*)
610	READ(inum,*)
611	READ(inum,*)
612	READ(inum,*) (bfoe(jj),jj=jpjed, jpjef)
613	CLOSE(inum)
614	END IF
615	DO jj=jpjed, jpjefm1
616	bfoe(jj)=bfoe(jj)*zcoef
617	END DO
618
619	END IF
620
621	IF( lpwestobc) THEN
622
623	IF( kt == nit000 .OR. kt <= kbsfstart ) THEN
624	OPEN(inum,file='obcwestbsf.dta')
625	READ(inum,*)
626	READ(inum,*)
627	READ(inum,*)
628	READ(inum,*)
629	READ(inum,*)
630	READ(inum,*) (bfow(jj),jj=jpjwd, jpjwf)
631	CLOSE(inum)
632	END IF
633	DO jj=jpjwd, jpjwfm1
634	bfow(jj)=bfow(jj)*zcoef
635	END DO
636
637	END IF
638
639	IF( lpsouthobc) THEN
640
641	IF( kt == nit000.OR.kt <= kbsfstart ) THEN
642	OPEN(inum,file='obcsouthbsf.dta')
643	READ(inum,*)
644	READ(inum,*)
645	READ(inum,*)
646	READ(inum,*)
647	READ(inum,*)
648	READ(inum,*) (bfos(jj),jj=jpisd, jpisf)
649	CLOSE(inum)
650	END IF
651	DO ji=jpisd, jpisfm1
652	bfos(ji)=bfos(ji)*zcoef
653	END DO
654
655	END IF
656
657	IF( lpnorthobc) THEN
658	IF( kt == nit000.OR.kt <= kbsfstart ) THEN
659	OPEN(inum,file='obcnorthbsf.dta')
660	READ(inum,*)
661	READ(inum,*)
662	READ(inum,*)
663	READ(inum,*)
664	READ(inum,*)
665	READ(inum,*) (bfon(jj),jj=jpind, jpinf)
666	CLOSE(inum)
667	END IF
668	DO ji=jpind, jpinfm1
669	bfon(ji)=bfon(ji)*zcoef
670	END DO
671
672	END IF
673	END SUBROUTINE obc_dta_psi
674
675	#else
676	!!-----------------------------------------------------------------------------
677	!! 'key_dynspg_rl not defined ( dynspg_flt or dynspg_exp or dynspg_ts )
678	!!-----------------------------------------------------------------------------
679	SUBROUTINE obc_dta_psi ( kt ) ! Empty routine
680	!! * Arguments
681	INTEGER,INTENT(in) :: kt
682	WRITE(,) 'obc_dta_psi: You should not have seen this print! error?', kt
683	STOP
684	END SUBROUTINE obc_dta_psi
685	# endif
686
687	#if defined key_dynspg_ts \|\| defined key_dynspg_exp
688	SUBROUTINE obc_dta_bt( kt, kbt )
689	!!---------------------------------------------------------------------------
690	!! * SUBROUTINE obc_dta *
691	!!
692	!! ** Purpose : time interpolation of barotropic data for time-splitting scheme
693	!! Data at the boundary must be in m2/s
694	!!
695	!! History :
696	!! 9.0 ! 05-11 (V. garnier) Original code
697	!!---------------------------------------------------------------------------
698	!! * Arguments
699	INTEGER, INTENT( in ) :: kt ! ocean time-step index
700	INTEGER, INTENT( in ) :: kbt ! barotropic ocean time-step index
701
702	!! * Local declarations
703	INTEGER :: ji, jj, jk, ii, ij ! dummy loop indices
704	INTEGER :: id_e, id_w, id_n, id_s, fid ! file identifiers
705	INTEGER :: itimo, iman, imois, i15
706	INTEGER :: itobcm, itobcp, itimom, itimop
707	REAL(wp) :: zxy
708	INTEGER :: isrel, ikt ! number of seconds since 1/1/1992
709	INTEGER :: iprint ! frequency for printouts.
710
711	!!---------------------------------------------------------------------------
712
713	! 1. First call: check time frames available in files.
714	! -------------------------------------------------------
715
716	IF( kt == nit000 ) THEN
717
718	! 1.1 Barotropic tangential velocities set to zero
719	! -------------------------------------------------
720	IF( lp_obc_east ) vbtfoe(:) = 0.e0
721	IF( lp_obc_west ) vbtfow(:) = 0.e0
722	IF( lp_obc_south ) ubtfos(:) = 0.e0
723	IF( lp_obc_north ) ubtfon(:) = 0.e0
724
725	! 1.2 Sea surface height and normal barotropic velocities set to zero
726	! or initial conditions if nobc_dta == 0
727	! --------------------------------------------------------------------
728
729	IF( lp_obc_east ) THEN
730	! initialisation to zero
731	sshedta(:,:) = 0.e0
732	ubtedta(:,:) = 0.e0
733	! ! ================== !
734	IF( nobc_dta == 0 ) THEN ! initial state used !
735	! ! ================== !
736	! Fills sedta, tedta, uedta (global arrays)
737	! Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
738	DO ji = nie0, nie1
739	DO jj = nje0p1, nje1m1
740	ij = jj -1 + njmpp
741	sshedta(ij,1) = sshn(ji+1,jj) * tmask(ji+1,jj,1)
742	END DO
743	END DO
744	ENDIF
745	ENDIF
746
747	IF( lp_obc_west) THEN
748	! initialisation to zero
749	sshwdta(:,:) = 0.e0
750	ubtwdta(:,:) = 0.e0
751	! ! ================== !
752	IF( nobc_dta == 0 ) THEN ! initial state used !
753	! ! ================== !
754	! Fills swdta, twdta, uwdta (global arrays)
755	! Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
756	DO ji = niw0, niw1
757	DO jj = njw0p1, njw1m1
758	ij = jj -1 + njmpp
759	sshwdta(ij,1) = sshn(ji,jj) * tmask(ji,jj,1)
760	END DO
761	END DO
762	ENDIF
763	ENDIF
764
765	IF( lp_obc_north) THEN
766	! initialisation to zero
767	sshndta(:,:) = 0.e0
768	vbtndta(:,:) = 0.e0
769	! ! ================== !
770	IF( nobc_dta == 0 ) THEN ! initial state used !
771	! ! ================== !
772	! Fills sndta, tndta, vndta (global arrays)
773	! Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
774	DO jj = njn0, njn1
775	DO ji = nin0p1, nin1m1
776	DO jk = 1, jpkm1
777	ii = ji -1 + nimpp
778	vbtndta(ii,1) = vbtndta(ii,1) + vndta(ii,jk,1)*fse3v(ji,jj,jk)
779	END DO
780	sshndta(ii,1) = sshn(ii,jj+1) * tmask(ji,jj+1,1)
781	END DO
782	END DO
783	ENDIF
784	ENDIF
785
786	IF( lp_obc_south) THEN
787	! initialisation to zero
788	ssdta(:,:,:) = 0.e0
789	tsdta(:,:,:) = 0.e0
790	vsdta(:,:,:) = 0.e0
791	sshsdta(:,:) = 0.e0
792	vbtsdta(:,:) = 0.e0
793	! ! ================== !
794	IF( nobc_dta == 0 ) THEN ! initial state used !
795	! ! ================== !
796	! Fills ssdta, tsdta, vsdta (global arrays)
797	! Remark: this works for njzoom = 1. Should the definition of ij include njzoom?
798	DO jj = njs0, njs1
799	DO ji = nis0p1, nis1m1
800	DO jk = 1, jpkm1
801	ii = ji -1 + nimpp
802	vbtsdta(ii,1) = vbtsdta(ii,1) + vsdta(ii,jk,1)*fse3v(ji,jj,jk)
803	END DO
804	sshsdta(ii,1) = sshn(ji,jj) * tmask(ii,jj,1)
805	END DO
806	END DO
807	ENDIF
808	ENDIF
809
810	ENDIF ! END IF kt == nit000
811
812	!!------------------------------------------------------------------------------------
813	! 2. Initialize the time we are at. Does this every time the routine is called,
814	! excepted when nobc_dta = 0
815	!
816	IF( nobc_dta == 0) THEN
817	itimo = 1
818	zxy = 0
819	ELSE
820	IF(ntobc == 1) THEN
821	itimo = 1
822	ELSE IF (ntobc == 12) THEN ! BC are monthly
823	! we assume we have climatology in that case
824	iman = 12
825	i15 = nday / 16
826	imois = nmonth + i15 - 1
827	IF( imois == 0 ) imois = iman
828	itimo = imois
829	ELSE
830	IF(lwp) WRITE(numout,*) 'data other than constant or monthly',kt
831	iman = ntobc
832	itimo = FLOOR( kt*rdt / tcobc(1))
833	isrel=kt*rdt
834	ENDIF
835	ENDIF
836
837	! 2. Read two records in the file if necessary
838	! ---------------------------------------------
839
840	IF( nobc_dta == 1 .AND. nlecto == 1 ) THEN
841
842	IF( lp_obc_east ) THEN
843	! ... Read datafile and set sea surface height and barotropic velocity
844	! ... initialise the sshedta, ubtedta arrays
845	sshedta(:,0) = sshedta(:,1)
846	ubtedta(:,0) = ubtedta(:,1)
847	CALL iom_open ( 'obceast_TS.nc', id_e )
848	CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,1), ktime=ntobc1 )
849	CALL iom_get ( id_e, jpdom_unknown, 'vossurfh', sshedta(:,2), ktime=ntobc2 )
850	IF( lk_dynspg_ts ) THEN
851	CALL iom_get (id_e, jpdom_unknown, 'vossurfh', sshedta(:,3), ktime=ntobc2+1 )
852	ENDIF
853	CALL iom_close ( id_e )
854	!
855	CALL iom_open ( 'obceast_U.nc', id_e )
856	CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,1), ktime=ntobc1 )
857	CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,2), ktime=ntobc2 )
858	IF( lk_dynspg_ts ) THEN
859	CALL iom_get ( id_e, jpdom_unknown, 'vozoubt', ubtedta(:,3), ktime=ntobc2+1 )
860	ENDIF
861	CALL iom_close ( id_e )
862	! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
863	IF( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
864	WRITE(numout,*)
865	WRITE(numout,*) ' Read East OBC barotropic data records ', ntobc1, ntobc2
866	iprint = (jpjef-jpjed+1)/20 +1
867	WRITE(numout,*)
868	WRITE(numout,*) ' Sea surface height record 1'
869	CALL prihre( sshedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout )
870	WRITE(numout,*)
871	WRITE(numout,*) ' Normal transport (m2/s) record 1',iprint
872	CALL prihre( ubtedta(:,1), jpjef-jpjed+1, 1, 1, jpjef-jpjed+1, iprint, 1, 1, -3, 1., numout )
873	ENDIF
874	ENDIF
875
876	IF( lp_obc_west ) THEN
877	! ... Read datafile and set temperature, salinity and normal velocity
878	! ... initialise the swdta, twdta, uwdta arrays
879	sshwdta(:,0) = sshwdta(:,1)
880	ubtwdta(:,0) = ubtwdta(:,1)
881	CALL iom_open ( 'obcwest_TS.nc', id_w )
882	CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,1), ktime=ntobc1 )
883	CALL iom_get ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,2), ktime=ntobc2 )
884	IF( lk_dynspg_ts ) THEN
885	CALL ( id_w, jpdom_unknown, 'vossurfh', sshwdta(:,3), ktime=ntobc2+1 )
886	ENDIF
887	CALL iom_close ( id_w )
888	!
889	CALL iom_open ( 'obcwest_U.nc', id_w )
890	CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,1), ktime=ntobc1 )
891	CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,2), ktime=ntobc2 )
892	IF( lk_dynspg_ts ) THEN
893	CALL iom_get ( id_w, jpdom_unknown, 'vozoubt', ubtwdta(:,3), ktime=ntobc2+1 )
894	ENDIF
895	CALL iom_close ( id_w )
896	! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
897	IF( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
898	WRITE(numout,*)
899	WRITE(numout,*) ' Read West OBC barotropic data records ', ntobc1, ntobc2
900	iprint = (jpjwf-jpjwd+1)/20 +1
901	WRITE(numout,*)
902	WRITE(numout,*) ' Sea surface height record 1 - printout surface level'
903	CALL prihre( sshwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout )
904	WRITE(numout,*)
905	WRITE(numout,*) ' Normal transport (m2/s) record 1'
906	CALL prihre( ubtwdta(:,1), jpjwf-jpjwd+1, 1, 1, jpjwf-jpjwd+1, iprint, 1, 1, -3, 1., numout )
907	ENDIF
908	ENDIF
909
910	IF( lp_obc_north) THEN
911	! ... Read datafile and set sea surface height and barotropic velocity
912	! ... initialise the sshndta, ubtndta arrays
913	sshndta(:,0) = sshndta(:,1)
914	vbtndta(:,0) = vbtndta(:,1)
915	CALL iom_open ( 'obcnorth_TS.nc', id_n )
916	CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,1), ktime=ntobc1 )
917	CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,2), ktime=ntobc2 )
918	IF( lk_dynspg_ts ) THEN
919	CALL iom_get (id_n, jpdom_unknown, 'vossurfh', sshndta(:,3), ktime=ntobc2+1 )
920	ENDIF
921	CALL iom_close ( id_n )
922
923	CALL iom_open ( 'obcnorth_V.nc', id_n )
924	CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,1), ktime=ntobc1 )
925	CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,2), ktime=ntobc2 )
926	IF( lk_dynspg_ts ) THEN
927	CALL iom_get (id_n, jpdom_unknown, 'vomevbt', vbtndta(:,3), ktime=ntobc2+1 )
928	ENDIF
929	CALL iom_close ( id_n )
930
931	! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
932	IF( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
933	WRITE(numout,*)
934	WRITE(numout,*) ' Read North OBC barotropic data records ', ntobc1, ntobc2
935	iprint = (jpinf-jpind+1)/20 +1
936	WRITE(numout,*)
937	WRITE(numout,*) ' Sea surface height record 1 - printout surface level'
938	CALL prihre( sshndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout )
939	WRITE(numout,*)
940	WRITE(numout,*) ' Normal transport (m2/s) record 1'
941	CALL prihre( vbtndta(:,1), jpinf-jpind+1, 1, 1, jpinf-jpind+1, iprint, 1, 1, -3, 1., numout )
942	ENDIF
943	ENDIF
944
945	IF( lp_obc_south) THEN
946	! ... Read datafile and set sea surface height and barotropic velocity
947	! ... initialise the sshsdta, ubtsdta arrays
948	sshsdta(:,0) = sshsdta(:,1)
949	vbtsdta(:,0) = vbtsdta(:,1)
950	CALL iom_open ( 'obcsouth_TS.nc', id_s )
951	CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,1), ktime=ntobc1 )
952	CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,2), ktime=ntobc2 )
953	IF( lk_dynspg_ts ) THEN
954	CALL iom_get ( id_s, jpdom_unknown, 'vossurfh', sshsdta(:,3), ktime=ntobc2+1 )
955	ENDIF
956	CALL iom_close ( id_s )
957
958	CALL iom_open ( 'obcsouth_V.nc', id_s )
959	CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,1), ktime=ntobc1 )
960	CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,2), ktime=ntobc2 )
961	IF( lk_dynspg_ts ) THEN
962	CALL iom_get ( id_s, jpdom_unknown, 'vomevbt', vbtsdta(:,3), ktime=ntobc2+1 )
963	ENDIF
964	CALL iom_close ( id_s )
965
966	! ... Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
967	IF( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
968	WRITE(numout,*)
969	WRITE(numout,*) ' Read South OBC barotropic data records ', ntobc1, ntobc2
970	iprint = (jpisf-jpisd+1)/20 +1
971	WRITE(numout,*)
972	WRITE(numout,*) ' Sea surface height record 1 - printout surface level'
973	CALL prihre( sshsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout )
974	WRITE(numout,*)
975	WRITE(numout,*) ' Normal transport (m2/s) record 1'
976	CALL prihre( vbtsdta(:,1), jpisf-jpisd+1, 1, 1, jpisf-jpisd+1, iprint, 1, 1, -3, 1., numout )
977	ENDIF
978	ENDIF
979
980	ENDIF ! end of the test on the condition to read or not the files
981
982	! 3. Call at every time step : Linear interpolation of BCs to current time step
983	! ----------------------------------------------------------------------
984
985	IF( lk_dynspg_ts ) THEN
986	isrel = (kt-1)rdt + kbtrdtbt
987
988	IF( nobc_dta == 1 ) THEN
989	isrel = (kt-1)rdt + kbtrdtbt
990	itimo = FLOOR( kt*rdt / (tcobc(2)-tcobc(1)) )
991	itimom = FLOOR( (kt-1)*rdt / (tcobc(2)-tcobc(1)) )
992	itimop = FLOOR( (kt+1)*rdt / (tcobc(2)-tcobc(1)) )
993	IF( itimom == itimo .AND. itimop == itimo ) THEN
994	itobcm = ntobc1
995	itobcp = ntobc2
996
997	ELSEIF ( itimom <= itimo .AND. itimop == itimo ) THEN
998	IF( FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN
999	itobcm = ntobc1-1
1000	itobcp = ntobc2-1
1001	ELSE
1002	itobcm = ntobc1
1003	itobcp = ntobc2
1004	ENDIF
1005
1006	ELSEIF ( itimom == itimo .AND. itimop >= itimo ) THEN
1007	IF( FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN
1008	itobcm = ntobc1
1009	itobcp = ntobc2
1010	ELSE
1011	itobcm = ntobc1+1
1012	itobcp = ntobc2+1
1013	ENDIF
1014
1015	ELSEIF ( itimom == itimo-1 .AND. itimop == itimo+1 ) THEN
1016	IF( FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimo ) THEN
1017	itobcm = ntobc1-1
1018	itobcp = ntobc2-1
1019	ELSEIF ( FLOOR( isrel / (tcobc(2)-tcobc(1)) ) < itimop ) THEN
1020	itobcm = ntobc1
1021	itobcp = ntobc2
1022	ELSEIF ( FLOOR( isrel / (tcobc(2)-tcobc(1)) ) == itimop ) THEN
1023	itobcm = ntobc1+1
1024	itobcp = ntobc2+2
1025	ELSE
1026	IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 1?'
1027	ENDIF
1028	ELSE
1029	IF(lwp) WRITE(numout, *) 'obc_dta_bt: You should not have seen this print! error 2?'
1030	ENDIF
1031
1032	ENDIF
1033
1034	ELSE IF( lk_dynspg_exp ) THEN
1035	isrel=kt*rdt
1036	itobcm = ntobc1
1037	itobcp = ntobc2
1038	ENDIF
1039
1040	IF( ntobc == 1 .OR. nobc_dta == 0 ) THEN
1041	zxy = 0.e0
1042	ELSE IF( ntobc == 12 ) THEN
1043	zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
1044	ELSE
1045	zxy = (tcobc(itobcm)-FLOAT(isrel)) / (tcobc(itobcm)-tcobc(itobcp))
1046	ENDIF
1047
1048	IF( lp_obc_east ) THEN ! fills sshfoe, ubtfoe (local to each processor)
1049	DO jj = nje0p1, nje1m1
1050	ij = jj -1 + njmpp
1051	sshfoe(jj) = ( zxy * sshedta(ij,2) + (1.-zxy) * sshedta(ij,1) ) * temsk(jj,1)
1052	ubtfoe(jj) = ( zxy * ubtedta(ij,2) + (1.-zxy) * ubtedta(ij,1) ) * uemsk(jj,1)
1053	END DO
1054	ENDIF
1055
1056	IF( lp_obc_west) THEN ! fills sshfow, ubtfow (local to each processor)
1057	DO jj = njw0p1, njw1m1
1058	ij = jj -1 + njmpp
1059	sshfow(jj) = ( zxy * sshwdta(ij,2) + (1.-zxy) * sshwdta(ij,1) ) * twmsk(jj,1)
1060	ubtfow(jj) = ( zxy * ubtwdta(ij,2) + (1.-zxy) * ubtwdta(ij,1) ) * uwmsk(jj,1)
1061	END DO
1062	ENDIF
1063
1064	IF( lp_obc_north) THEN ! fills sshfon, vbtfon (local to each processor)
1065	DO ji = nin0p1, nin1m1
1066	ii = ji -1 + nimpp
1067	sshfon(ji) = ( zxy * sshndta(ii,2) + (1.-zxy) * sshndta(ii,1) ) * tnmsk(ji,1)
1068	vbtfon(ji) = ( zxy * vbtndta(ii,2) + (1.-zxy) * vbtndta(ii,1) ) * vnmsk(ji,1)
1069	END DO
1070	ENDIF
1071
1072	IF( lp_obc_south) THEN ! fills sshfos, vbtfos (local to each processor)
1073	DO ji = nis0p1, nis1m1
1074	ii = ji -1 + nimpp
1075	sshfos(ji) = ( zxy * sshsdta(ii,2) + (1.-zxy) * sshsdta(ii,1) ) * tsmsk(ji,1)
1076	vbtfos(ji) = ( zxy * vbtsdta(ii,2) + (1.-zxy) * vbtsdta(ii,1) ) * vsmsk(ji,1)
1077	END DO
1078	ENDIF
1079
1080	END SUBROUTINE obc_dta_bt
1081
1082	#else
1083	!!-----------------------------------------------------------------------------
1084	!! Default option
1085	!!-----------------------------------------------------------------------------
1086	SUBROUTINE obc_dta_bt ( kt, kbt ) ! Empty routine
1087	!! * Arguments
1088	INTEGER,INTENT(in) :: kt
1089	INTEGER, INTENT( in ) :: kbt ! barotropic ocean time-step index
1090	WRITE(,) 'obc_dta_bt: You should not have seen this print! error?', kt
1091	WRITE(,) 'obc_dta_bt: You should not have seen this print! error?', kbt
1092	END SUBROUTINE obc_dta_bt
1093	#endif
1094
1095
1096	!!==============================================================================
1097	SUBROUTINE obc_read (kt, nt_x, ntobc_x, iyy, imm)
1098	!!-------------------------------------------------------------------------
1099	!! * ROUTINE obc_read *
1100	!!
1101	!! ** Purpose : Read the boundary data in files identified by iyy and imm
1102	!! According to the validated open boundaries, return the
1103	!! following arrays :
1104	!! sedta, tedta : East OBC salinity and temperature
1105	!! uedta, vedta : " " u and v velocity component
1106	!!
1107	!! swdta, twdta : West OBC salinity and temperature
1108	!! uwdta, vwdta : " " u and v velocity component
1109	!!
1110	!! sndta, tndta : North OBC salinity and temperature
1111	!! undta, vndta : " " u and v velocity component
1112	!!
1113	!! ssdta, tsdta : South OBC salinity and temperature
1114	!! usdta, vsdta : " " u and v velocity component
1115	!!
1116	!! ** Method : These fields are read in the record ntobc_x of the files.
1117	!! The number of records is already known. If ntobc_x is greater
1118	!! than the number of record, this routine will look for next file,
1119	!! updating the indices (case of inter-annual obcs) or loop at the
1120	!! begining in case of climatological file (ln_obc_clim = true ).
1121	!! -------------------------------------------------------------------------
1122	!! History: ! 2005 ( P. Mathiot, C. Langlais ) Original code
1123	!! ! 2008 ( J,M, Molines ) Use IOM and cleaning
1124	!!--------------------------------------------------------------------------
1125
1126	! * Arguments
1127	INTEGER, INTENT( in ) :: kt, nt_x
1128	INTEGER, INTENT( inout ) :: ntobc_x , iyy, imm ! yes ! inout !
1129
1130	! * Local variables
1131	CHARACTER (len=40) :: & ! file names
1132	cl_obc_eTS , cl_obc_eU, cl_obc_eV,&
1133	cl_obc_wTS , cl_obc_wU, cl_obc_wV,&
1134	cl_obc_nTS , cl_obc_nU, cl_obc_nV,&
1135	cl_obc_sTS , cl_obc_sU, cl_obc_sV
1136
1137	INTEGER :: ikprint
1138	REAL(wp) :: zmin, zmax ! control of boundary values
1139
1140	!IOM stuff
1141	INTEGER :: id_e, id_w, id_n, id_s, ji, jj
1142	INTEGER, DIMENSION(2) :: istart, icount
1143
1144	!--------------------------------------------------------------------------
1145	IF ( ntobc_x > itobc ) THEN
1146	IF (ln_obc_clim) THEN ! just loop on the same file
1147	ntobc_x = 1
1148	ELSE
1149	! need to change file : it is always for an 'after' data
1150	IF ( cffile == 'annual' ) THEN ! go to next year file
1151	iyy = iyy + 1
1152	ELSE IF ( cffile =='monthly' ) THEN ! go to next month file
1153	imm = imm + 1
1154	IF ( imm == 13 ) THEN
1155	imm = 1 ; iyy = iyy + 1
1156	ENDIF
1157	ELSE
1158	ctmp1='obcread : this type of obc file is not supported :( '
1159	ctmp2=TRIM(cffile)
1160	CALL ctl_stop (ctmp1, ctmp2)
1161	! cffile should be either annual or monthly ...
1162	ENDIF
1163	! as the file is changed, need to update itobc etc ...
1164	CALL obc_dta_chktime (iyy,imm)
1165	ntobc_x = nrecbef() + 1 ! remember : this case occur for an after data
1166	ENDIF
1167	ENDIF
1168
1169	IF ( lp_obc_east ) THEN
1170	! ... Read datafile and set temperature, salinity and normal velocity
1171	! ... initialise the sedta, tedta, uedta arrays
1172	WRITE(cl_obc_eTS ,'("obc_east_TS_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1173	WRITE(cl_obc_eU ,'("obc_east_U_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1174	WRITE(cl_obc_eV ,'("obc_east_V_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1175	! JMM this may change depending on the obc data format ...
1176	istart(:)=(/nje0+njmpp-1,1/) ; icount(:)=(/nje1-nje0 +1,jpk/)
1177	IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_eTS)
1178	IF (nje1 >= nje0 ) THEN
1179	CALL iom_open ( cl_obc_eTS , id_e )
1180	CALL iom_get ( id_e, jpdom_unknown, 'votemper', tedta(nje0:nje1,:,nt_x), &
1181	& ktime=ntobc_x , kstart=istart, kcount= icount )
1182	CALL iom_get ( id_e, jpdom_unknown, 'vosaline', sedta(nje0:nje1,:,nt_x), &
1183	& ktime=ntobc_x , kstart=istart, kcount= icount )
1184	CALL iom_close (id_e)
1185	!
1186	CALL iom_open ( cl_obc_eU , id_e )
1187	CALL iom_get ( id_e, jpdom_unknown, 'vozocrtx', uedta(nje0:nje1,:,nt_x), &
1188	& ktime=ntobc_x , kstart=istart, kcount= icount )
1189	CALL iom_close ( id_e )
1190	!
1191	CALL iom_open ( cl_obc_eV , id_e )
1192	CALL iom_get ( id_e, jpdom_unknown, 'vomecrty', vedta(nje0:nje1,:,nt_x), &
1193	& ktime=ntobc_x , kstart=istart, kcount= icount )
1194	CALL iom_close ( id_e )
1195
1196	! mask the boundary values
1197	tedta(:,:,nt_x) = tedta(:,:,nt_x)temsk(:,:) ; sedta(:,:,nt_x) = sedta(:,:,nt_x)temsk(:,:)
1198	uedta(:,:,nt_x) = uedta(:,:,nt_x)uemsk(:,:) ; vedta(:,:,nt_x) = vedta(:,:,nt_x)vemsk(:,:)
1199
1200	! check any outliers
1201	zmin=MINVAL( sedta(:,:,nt_x), mask=ltemsk ) ; zmax=MAXVAL(sedta(:,:,nt_x), mask=ltemsk)
1202	IF ( zmin < 5 .OR. zmax > 50) THEN
1203	CALL ctl_stop('Error in sedta',' routine obcdta')
1204	ENDIF
1205	zmin=MINVAL( tedta(:,:,nt_x), mask=ltemsk ) ; zmax=MAXVAL(tedta(:,:,nt_x), mask=ltemsk)
1206	IF ( zmin < -10. .OR. zmax > 40) THEN
1207	CALL ctl_stop('Error in tedta',' routine obcdta')
1208	ENDIF
1209	zmin=MINVAL( uedta(:,:,nt_x), mask=luemsk ) ; zmax=MAXVAL(uedta(:,:,nt_x), mask=luemsk)
1210	IF ( zmin < -5. .OR. zmax > 5.) THEN
1211	CALL ctl_stop('Error in uedta',' routine obcdta')
1212	ENDIF
1213	zmin=MINVAL( vedta(:,:,nt_x), mask=lvemsk ) ; zmax=MAXVAL(vedta(:,:,nt_x), mask=lvemsk)
1214	IF ( zmin < -5. .OR. zmax > 5.) THEN
1215	CALL ctl_stop('Error in vedta',' routine obcdta')
1216	ENDIF
1217
1218	! Usually printout is done only once at kt = nit000, unless nprint (namelist) > 1
1219	IF ( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
1220	WRITE(numout,*)
1221	WRITE(numout,*) ' Read East OBC data records ', ntobc_x
1222	ikprint = jpj/20 +1
1223	WRITE(numout,*) ' Temperature record 1 - printout every 3 level'
1224	CALL prihre( tedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1225	WRITE(numout,*)
1226	WRITE(numout,*) ' Salinity record 1 - printout every 3 level'
1227	CALL prihre( sedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1228	WRITE(numout,*)
1229	WRITE(numout,*) ' Normal velocity U record 1 - printout every 3 level'
1230	CALL prihre( uedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1231	WRITE(numout,*)
1232	WRITE(numout,*) ' Tangential velocity V record 1 - printout every 3 level'
1233	CALL prihre( vedta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1234	ENDIF
1235	ENDIF
1236	ENDIF
1237	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1238	IF ( lp_obc_west ) THEN
1239	! ... Read datafile and set temperature, salinity and normal velocity
1240	! ... initialise the swdta, twdta, uwdta arrays
1241	WRITE(cl_obc_wTS ,'("obc_west_TS_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1242	WRITE(cl_obc_wU ,'("obc_west_U_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1243	WRITE(cl_obc_wV ,'("obc_west_V_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1244	istart(:)=(/njw0+njmpp-1,1/) ; icount(:)=(/njw1-njw0 +1,jpk/)
1245	IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_wTS)
1246
1247	IF ( njw1 >= njw0 ) THEN
1248	CALL iom_open ( cl_obc_wTS , id_w )
1249	CALL iom_get ( id_w, jpdom_unknown, 'votemper', twdta(njw0:njw1,:,nt_x), &
1250	& ktime=ntobc_x , kstart=istart, kcount= icount )
1251
1252	CALL iom_get ( id_w, jpdom_unknown, 'vosaline', swdta(njw0:njw1,:,nt_x), &
1253	& ktime=ntobc_x , kstart=istart, kcount= icount)
1254	CALL iom_close (id_w)
1255	!
1256	CALL iom_open ( cl_obc_wU , id_w )
1257	CALL iom_get ( id_w, jpdom_unknown, 'vozocrtx', uwdta(njw0:njw1,:,nt_x),&
1258	& ktime=ntobc_x , kstart=istart, kcount= icount )
1259	CALL iom_close ( id_w )
1260	!
1261	CALL iom_open ( cl_obc_wV , id_w )
1262	CALL iom_get ( id_w, jpdom_unknown, 'vomecrty', vwdta(njw0:njw1,:,nt_x), &
1263	& ktime=ntobc_x , kstart=istart, kcount= icount )
1264	CALL iom_close ( id_w )
1265
1266	! mask the boundary values
1267	twdta(:,:,nt_x) = twdta(:,:,nt_x)twmsk(:,:) ; swdta(:,:,nt_x) = swdta(:,:,nt_x)twmsk(:,:)
1268	uwdta(:,:,nt_x) = uwdta(:,:,nt_x)uwmsk(:,:) ; vwdta(:,:,nt_x) = vwdta(:,:,nt_x)vwmsk(:,:)
1269
1270	! check any outliers
1271	zmin=MINVAL( swdta(:,:,nt_x), mask=ltwmsk ) ; zmax=MAXVAL(swdta(:,:,nt_x), mask=ltwmsk)
1272	IF ( zmin < 5 .OR. zmax > 50) THEN
1273	CALL ctl_stop('Error in swdta',' routine obcdta')
1274	ENDIF
1275	zmin=MINVAL( twdta(:,:,nt_x), mask=ltwmsk ) ; zmax=MAXVAL(twdta(:,:,nt_x), mask=ltwmsk)
1276	IF ( zmin < -10. .OR. zmax > 40) THEN
1277	CALL ctl_stop('Error in twdta',' routine obcdta')
1278	ENDIF
1279	zmin=MINVAL( uwdta(:,:,nt_x), mask=luwmsk ) ; zmax=MAXVAL(uwdta(:,:,nt_x), mask=luwmsk)
1280	IF ( zmin < -5. .OR. zmax > 5.) THEN
1281	CALL ctl_stop('Error in uwdta',' routine obcdta')
1282	ENDIF
1283	zmin=MINVAL( vwdta(:,:,nt_x), mask=lvwmsk ) ; zmax=MAXVAL(vwdta(:,:,nt_x), mask=lvwmsk)
1284	IF ( zmin < -5. .OR. zmax > 5.) THEN
1285	CALL ctl_stop('Error in vwdta',' routine obcdta')
1286	ENDIF
1287
1288
1289	IF ( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
1290	WRITE(numout,*)
1291	WRITE(numout,*) ' Read West OBC data records ', ntobc_x
1292	ikprint = jpj/20 +1
1293	WRITE(numout,*) ' Temperature record 1 - printout every 3 level'
1294	CALL prihre( twdta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1295	WRITE(numout,*)
1296	WRITE(numout,*) ' Salinity record 1 - printout every 3 level'
1297	CALL prihre( swdta(:,:,nt_x),jpj,jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1298	WRITE(numout,*)
1299	WRITE(numout,*) ' Normal velocity U record 1 - printout every 3 level'
1300	CALL prihre( uwdta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1301	WRITE(numout,*)
1302	WRITE(numout,*) ' Tangential velocity V record 1 - printout every 3 level'
1303	CALL prihre( vwdta(:,:,nt_x), jpj, jpk, 1, jpj, ikprint, jpk, 1, -3, 1., numout )
1304	ENDIF
1305	END IF
1306	ENDIF
1307	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1308	IF( lp_obc_north) THEN
1309	WRITE(cl_obc_nTS ,'("obc_north_TS_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1310	WRITE(cl_obc_nV ,'("obc_north_V_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1311	WRITE(cl_obc_nU ,'("obc_north_U_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1312	istart(:)=(/nin0+nimpp-1,1/) ; icount(:)=(/nin1-nin0 +1,jpk/)
1313	IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_nTS)
1314	IF ( nin1 >= nin0 ) THEN
1315	CALL iom_open ( cl_obc_nTS , id_n )
1316	CALL iom_get ( id_n, jpdom_unknown, 'votemper', tndta(nin0:nin1,:,nt_x), &
1317	& ktime=ntobc_x , kstart=istart, kcount= icount )
1318	CALL iom_get ( id_n, jpdom_unknown, 'vosaline', sndta(nin0:nin1,:,nt_x), &
1319	& ktime=ntobc_x , kstart=istart, kcount= icount )
1320	CALL iom_close (id_n)
1321	!
1322	CALL iom_open ( cl_obc_nU , id_n )
1323	CALL iom_get ( id_n, jpdom_unknown, 'vozocrtx', undta(nin0:nin1,:,nt_x), &
1324	& ktime=ntobc_x , kstart=istart, kcount= icount )
1325	CALL iom_close ( id_n )
1326	!
1327	CALL iom_open ( cl_obc_nV , id_n )
1328	CALL iom_get ( id_n, jpdom_unknown, 'vomecrty', vndta(nin0:nin1,:,nt_x), &
1329	& ktime=ntobc_x , kstart=istart, kcount= icount )
1330	CALL iom_close ( id_n )
1331
1332	! mask the boundary values
1333	tndta(:,:,nt_x) = tndta(:,:,nt_x)tnmsk(:,:) ; sndta(:,:,nt_x) = sndta(:,:,nt_x)tnmsk(:,:)
1334	undta(:,:,nt_x) = undta(:,:,nt_x)unmsk(:,:) ; vndta(:,:,nt_x) = vndta(:,:,nt_x)vnmsk(:,:)
1335
1336	! check any outliers
1337	zmin=MINVAL( sndta(:,:,nt_x), mask=ltnmsk ) ; zmax=MAXVAL(sndta(:,:,nt_x), mask=ltnmsk)
1338	IF ( zmin < 5 .OR. zmax > 50) THEN
1339	CALL ctl_stop('Error in sndta',' routine obcdta')
1340	ENDIF
1341	zmin=MINVAL( tndta(:,:,nt_x), mask=ltnmsk ) ; zmax=MAXVAL(tndta(:,:,nt_x), mask=ltnmsk)
1342	IF ( zmin < -10. .OR. zmax > 40) THEN
1343	CALL ctl_stop('Error in tndta',' routine obcdta')
1344	ENDIF
1345	zmin=MINVAL( undta(:,:,nt_x), mask=lunmsk ) ; zmax=MAXVAL(undta(:,:,nt_x), mask=lunmsk)
1346	IF ( zmin < -5. .OR. zmax > 5.) THEN
1347	CALL ctl_stop('Error in undta',' routine obcdta')
1348	ENDIF
1349	zmin=MINVAL( vndta(:,:,nt_x), mask=lvnmsk ) ; zmax=MAXVAL(vndta(:,:,nt_x), mask=lvnmsk)
1350	IF ( zmin < -5. .OR. zmax > 5.) THEN
1351	CALL ctl_stop('Error in vndta',' routine obcdta')
1352	ENDIF
1353
1354	IF ( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
1355	WRITE(numout,*)
1356	WRITE(numout,*) ' Read North OBC data records ', ntobc_x
1357	ikprint = jpi/20 +1
1358	WRITE(numout,*) ' Temperature record 1 - printout every 3 level'
1359	CALL prihre( tndta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1360	WRITE(numout,*)
1361	WRITE(numout,*) ' Salinity record 1 - printout every 3 level'
1362	CALL prihre( sndta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1363	WRITE(numout,*)
1364	WRITE(numout,*) ' Normal velocity V record 1 - printout every 3 level'
1365	CALL prihre( vndta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1366	WRITE(numout,*)
1367	WRITE(numout,*) ' Tangential velocity U record 1 - printout every 3 level'
1368	CALL prihre( undta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1369	ENDIF
1370	ENDIF
1371	ENDIF
1372	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1373	IF( lp_obc_south) THEN
1374	WRITE(cl_obc_sTS ,'("obc_south_TS_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1375	WRITE(cl_obc_sV ,'("obc_south_V_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1376	WRITE(cl_obc_sU ,'("obc_south_U_y" ,i4.4,"m",i2.2,".nc")' ) iyy,imm
1377	istart(:)=(/nis0+nimpp-1,1/) ; icount(:)=(/nis1-nis0 +1,jpk/)
1378	IF (lwp) WRITE(numout,*) 'read data in :', TRIM(cl_obc_sTS)
1379	IF ( nis1 >= nis0 ) THEN
1380	CALL iom_open ( cl_obc_sTS , id_s )
1381	CALL iom_get ( id_s, jpdom_unknown, 'votemper', tsdta(nis0:nis1,:,nt_x), &
1382	& ktime=ntobc_x , kstart=istart, kcount= icount )
1383	CALL iom_get ( id_s, jpdom_unknown, 'vosaline', ssdta(nis0:nis1,:,nt_x), &
1384	& ktime=ntobc_x , kstart=istart, kcount= icount )
1385	CALL iom_close (id_s)
1386	!
1387	CALL iom_open ( cl_obc_sU , id_s )
1388	CALL iom_get ( id_s, jpdom_unknown, 'vozocrtx', usdta(nis0:nis1,:,nt_x), &
1389	& ktime=ntobc_x , kstart=istart, kcount= icount )
1390	CALL iom_close ( id_s )
1391	!
1392	CALL iom_open ( cl_obc_sV , id_s )
1393	CALL iom_get ( id_s, jpdom_unknown, 'vomecrty', vsdta(nis0:nis1,:,nt_x), &
1394	& ktime=ntobc_x , kstart=istart, kcount= icount )
1395	CALL iom_close ( id_s )
1396
1397	! mask the boundary values
1398	tsdta(:,:,nt_x) = tsdta(:,:,nt_x)tsmsk(:,:) ; ssdta(:,:,nt_x) = ssdta(:,:,nt_x)tsmsk(:,:)
1399	usdta(:,:,nt_x) = usdta(:,:,nt_x)usmsk(:,:) ; vsdta(:,:,nt_x) = vsdta(:,:,nt_x)vsmsk(:,:)
1400
1401	! check any outliers
1402	zmin=MINVAL( ssdta(:,:,nt_x), mask=ltsmsk ) ; zmax=MAXVAL(ssdta(:,:,nt_x), mask=ltsmsk)
1403	IF ( zmin < 5 .OR. zmax > 50) THEN
1404	CALL ctl_stop('Error in ssdta',' routine obcdta')
1405	ENDIF
1406	zmin=MINVAL( tsdta(:,:,nt_x), mask=ltsmsk ) ; zmax=MAXVAL(tsdta(:,:,nt_x), mask=ltsmsk)
1407	IF ( zmin < -10. .OR. zmax > 40) THEN
1408	CALL ctl_stop('Error in tsdta',' routine obcdta')
1409	ENDIF
1410	zmin=MINVAL( usdta(:,:,nt_x), mask=lusmsk ) ; zmax=MAXVAL(usdta(:,:,nt_x), mask=lusmsk)
1411	IF ( zmin < -5. .OR. zmax > 5.) THEN
1412	CALL ctl_stop('Error in usdta',' routine obcdta')
1413	ENDIF
1414	zmin=MINVAL( vsdta(:,:,nt_x), mask=lvsmsk ) ; zmax=MAXVAL(vsdta(:,:,nt_x), mask=lvsmsk)
1415	IF ( zmin < -5. .OR. zmax > 5.) THEN
1416	CALL ctl_stop('Error in vsdta',' routine obcdta')
1417	ENDIF
1418
1419	IF ( lwp .AND. ( kt == nit000 .OR. nprint /= 0 ) ) THEN
1420	WRITE(numout,*)
1421	WRITE(numout,*) ' Read South OBC data records ', ntobc_x
1422	ikprint = jpi/20 +1
1423	WRITE(numout,*) ' Temperature record 1 - printout every 3 level'
1424	CALL prihre( tsdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1425	WRITE(numout,*)
1426	WRITE(numout,*) ' Salinity record 1 - printout every 3 level'
1427	CALL prihre( ssdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1428	WRITE(numout,*)
1429	WRITE(numout,*) ' Normal velocity V record 1 - printout every 3 level'
1430	CALL prihre( vsdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1431	WRITE(numout,*)
1432	WRITE(numout,*) ' Tangential velocity U record 1 - printout every 3 level'
1433	CALL prihre( usdta(:,:,nt_x), jpi, jpk, 1, jpi, ikprint, jpk, 1, -3, 1., numout )
1434	ENDIF
1435	ENDIF
1436	ENDIF
1437	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1438	END SUBROUTINE obc_read
1439
1440	INTEGER FUNCTION nrecbef()
1441	!!-----------------------------------------------------------------------
1442	!! * FUNCTION nrecbef *
1443	!!
1444	!! Purpose : - provide the before record number in files, with respect to zjcnes
1445	!!
1446	!! History : 2008-04 : ( J.M. Molines ) Original code
1447	!!-----------------------------------------------------------------------
1448
1449	INTEGER :: it , idum
1450
1451	idum = itobc
1452	DO it =1, itobc
1453	IF ( ztcobc(it) > zjcnes ) THEN ; idum = it - 1 ; EXIT ; ENDIF
1454	ENDDO
1455	! idum can be 0 (climato, before first record)
1456	IF ( idum == 0 ) THEN
1457	IF ( ln_obc_clim ) THEN
1458	idum = itobc
1459	ELSE
1460	ctmp1='obc_dta: find ntobc == 0 for non climatological file '
1461	ctmp2='consider adding a first record in your data file '
1462	CALL ctl_stop(ctmp1, ctmp2)
1463	ENDIF
1464	ENDIF
1465	! idum can be itobc ( zjcnes > ztcobc (itobc) )
1466	! This is not a problem ...
1467	nrecbef = idum
1468
1469	END FUNCTION nrecbef
1470
1471	#else
1472	!!------------------------------------------------------------------------------
1473	!! default option: Dummy module NO Open Boundary Conditions
1474	!!------------------------------------------------------------------------------
1475	CONTAINS
1476	SUBROUTINE obc_dta( kt ) ! Dummy routine
1477	INTEGER, INTENT (in) :: kt
1478	WRITE(,) 'obc_dta: You should not have seen this print! error?', kt
1479	END SUBROUTINE obc_dta
1480	#endif
1481	END MODULE obcdta

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

Download in other formats: