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

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

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

Download in other formats: