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

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

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

Download in other formats: