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

Last change on this file since 642 was 623, checked in by opalod, 17 years ago
nemo_v2_bugfix_024:RB: Use of ctlopn for all files (except for dimg and coupled part)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 52.3 KB

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

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

Download in other formats: