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

Last change on this file since 465 was 465, checked in by opalod, 18 years ago
nemo_v1_update_056:RB: updatye OBC with new coordinate
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 56.3 KB

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

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

Download in other formats: