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

Last change on this file since 389 was 389, checked in by opalod, 18 years ago

RB:nemo_v1_update_038: first integration of Agrif :

configuration parameters are just integer when agrif is used
add call to agrif routines with key_agrif

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats: