Context Navigation

Back to Ticket #377

Ticket #377: bdydta.F90

File bdydta.F90, 43.2 KB (added by ed.blockley, 15 years ago)
time offset calculation removed from bdydta.F90

Line
1	MODULE bdydta
2	!!======================================================================
3	!! * MODULE bdydta *
4	!! Open boundary data : read the data for the unstructured open boundaries.
5	!!======================================================================
6	!! History : 1.0 ! 2005-01 (J. Chanut, A. Sellar) Original code
7	!! - ! 2007-01 (D. Storkey) Update to use IOM module
8	!! - ! 2007-07 (D. Storkey) add bdy_dta_bt
9	!! 3.0 ! 2008-04 (NEMO team) add in the reference version
10	!!----------------------------------------------------------------------
11	#if defined key_bdy
12	!!----------------------------------------------------------------------
13	!! 'key_bdy' Unstructured Open Boundary Conditions
14	!!----------------------------------------------------------------------
15	!! bdy_dta : read u, v, t, s data along open boundaries
16	!! bdy_dta_bt : read depth-mean velocities and elevation along open
17	!! boundaries
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers
20	USE dom_oce ! ocean space and time domain
21	USE daymod ! calendar
22	USE phycst ! physical constants
23	USE bdy_oce ! ocean open boundary conditions
24	USE bdytides ! tidal forcing at boundaries
25	USE iom
26	USE ioipsl
27	USE in_out_manager ! I/O logical units
28
29	IMPLICIT NONE
30	PRIVATE
31
32	PUBLIC bdy_dta ! routines called by step.F90
33	PUBLIC bdy_dta_bt
34
35	INTEGER :: numbdyt, numbdyu, numbdyv !: logical units for T-, U-, & V-points data file, resp.
36	INTEGER :: ntimes_bdy !: exact number of time dumps in data files
37	INTEGER :: nbdy_b, nbdy_a !: record of bdy data file for before and after model time step
38	INTEGER :: numbdyt_bt, numbdyu_bt, numbdyv_bt !: logical unit for T-, U- & V-points data file, resp.
39	INTEGER :: ntimes_bdy_bt !: exact number of time dumps in data files
40	INTEGER :: nbdy_b_bt, nbdy_a_bt !: record of bdy data file for before and after model time step
41
42	INTEGER, DIMENSION (jpbtime) :: istep, istep_bt !: time array in seconds in each data file
43
44	REAL(wp) :: zoffset !: time offset between time origin in file & start time of model run
45
46	REAL(wp), DIMENSION(jpbdim,jpk,2) :: tbdydta, sbdydta !: time interpolated values of T and S bdy data
47	REAL(wp), DIMENSION(jpbdim,jpk,2) :: ubdydta, vbdydta !: time interpolated values of U and V bdy data
48	REAL(wp), DIMENSION(jpbdim,2) :: ubtbdydta, vbtbdydta !: Arrays used for time interpolation of bdy data
49	REAL(wp), DIMENSION(jpbdim,2) :: sshbdydta !: bdy data of ssh
50
51	!!----------------------------------------------------------------------
52	!! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
53	!! $Id: bdydta.F90 1241 2009-01-05 16:14:49Z rblod $
54	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
55	!!----------------------------------------------------------------------
56
57	CONTAINS
58
59	SUBROUTINE bdy_dta( kt )
60	!!----------------------------------------------------------------------
61	!! * SUBROUTINE bdy_dta *
62	!!
63	!! ** Purpose : Read unstructured boundary data for FRS condition.
64	!!
65	!! ** Method : At the first timestep, read in boundary data for two
66	!! times from the file and time-interpolate. At other
67	!! timesteps, check to see if we need another time from
68	!! the file. If so read it in. Time interpolate.
69	!!----------------------------------------------------------------------
70	INTEGER, INTENT( in ) :: kt ! ocean time-step index (for timesplitting option, otherwise zero)
71	!!
72	CHARACTER(LEN=80), DIMENSION(3) :: clfile ! names of input files
73	LOGICAL :: lect ! flag for reading
74	INTEGER :: it, ib, ik, igrd ! dummy loop indices
75	INTEGER :: igrd_start, igrd_end ! start and end of loops on igrd
76	INTEGER :: idvar ! netcdf var ID
77	INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing
78	INTEGER :: ntimes_bdyt, ntimes_bdyu, ntimes_bdyv
79	INTEGER :: itimer, totime
80	INTEGER :: ii, ij ! array addresses
81	INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read)
82	INTEGER :: iyear0, imonth0, iday0
83	INTEGER :: ihours0, iminutes0, isec0
84	INTEGER :: iyear, imonth, iday, isecs
85	INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files
86	REAL(wp) :: zxy, zoffsett
87	REAL(wp) :: zoffsetu, zoffsetv
88	REAL(wp) :: dayjul0, zdayjulini
89	REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files
90	REAL(wp), DIMENSION(jpbdta,1,jpk) :: zdta ! temporary array for data fields
91	!!---------------------------------------------------------------------------
92
93	IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs ) THEN ! If these are both false then this routine
94	! does nothing.
95
96	! -------------------- !
97	! Initialization !
98	! -------------------- !
99
100	lect = .false. ! If true, read a time record
101
102	! Some time variables for monthly climatological forcing:
103	! *******************************************************
104	!!gm here use directely daymod variables
105
106	iman = INT( raamo ) ! Number of months in a year
107
108	i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
109	! i15=0 if the current day is in the first half of the month, else i15=1
110
111	imois = nmonth + i15 - 1 ! imois is the first month record
112	IF( imois == 0 ) imois = iman
113
114	! Time variable for non-climatological forcing:
115	! *********************************************
116	itimer = (kt-nit000+1)*rdt ! current time in seconds for interpolation
117
118
119	! !-------------------!
120	IF( kt == nit000 ) THEN ! First call only !
121	! !-------------------!
122	istep(:) = 0
123	nbdy_b = 0
124	nbdy_a = 0
125
126	! Get time information from bdy data file
127	! ***************************************
128
129	IF(lwp) WRITE(numout,*)
130	IF(lwp) WRITE(numout,*) 'bdy_dta : Initialize unstructured boundary data'
131	IF(lwp) WRITE(numout,*) '~~~~~~~'
132
133	IF ( nbdy_dta == 0 ) THEN
134	!
135	IF(lwp) WRITE(numout,*) ' Bdy data are taken from initial conditions'
136	!
137	ELSEIF (nbdy_dta == 1) THEN
138	!
139	IF(lwp) WRITE(numout,*) ' Bdy data are read in netcdf files'
140	!
141	totime = ( nitend - nit000 + 1 ) * rdt ! Total time of the run to verify that all the
142	! ! necessary time dumps in file are included
143	!
144	clfile(1) = filbdy_data_T
145	clfile(2) = filbdy_data_U
146	clfile(3) = filbdy_data_V
147	!
148	! how many files are we to read in?
149	igrd_start = 1
150	igrd_end = 3
151	IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN
152	! No T-grid file.
153	igrd_start = 2
154	ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN
155	! No U-grid or V-grid file.
156	igrd_end = 1
157	ENDIF
158
159	DO igrd = igrd_start, igrd_end ! loop over T, U & V grid !
160	! !---------------------------!
161	CALL iom_open( clfile(igrd), inum )
162	CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, offset=zoffset, &
163	zyear=nyear, zmonth=nmonth, zday=nday, zsec_day=0.0 )
164
165	!! TO BE DONE... Check consistency between calendar from file
166	!! (available optionally from iom_gettime) and calendar in model
167	!! when calendar in model available outside of IOIPSL.
168
169	IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy
170	IF(lwp) WRITE(numout,*) 'offset: ',zoffset
171	IF(lwp) WRITE(numout,*) 'totime: ',totime
172	IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr
173
174	! Check that there are not too many times in the file.
175	IF( ntimes_bdy > jpbtime ) THEN
176	WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy
177	CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 )
178	ENDIF
179
180	! Check that time array increases:
181
182	it = 1
183	DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 )
184	it = it + 1
185	END DO
186
187	IF( it.NE.ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN
188	WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
189	CALL ctl_stop( 'Time array in unstructured boundary data files', &
190	& 'does not continuously increase.' , ctmp1 )
191	ENDIF
192	!
193	! Check that times in file span model run time:
194	IF( zstepr(1) + zoffset > 0 ) THEN
195	WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
196	CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 )
197	END IF
198	IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN
199	WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
200	CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 )
201	END IF
202	!
203	IF ( igrd == 1 ) THEN
204	ntimes_bdyt = ntimes_bdy
205	zoffsett = zoffset
206	istept(:) = INT( zstepr(:) + zoffset )
207	ELSEIF(igrd == 2 ) THEN
208	ntimes_bdyu = ntimes_bdy
209	zoffsetu = zoffset
210	istepu(:) = INT( zstepr(:) + zoffset )
211	ELSEIF(igrd == 3 ) THEN
212	ntimes_bdyv = ntimes_bdy
213	zoffsetv = zoffset
214	istepv(:) = INT( zstepr(:) + zoffset )
215	ENDIF
216	!
217	END DO ! end loop over T, U & V grid
218
219	IF (igrd_start == 1 .and. igrd_end == 3) THEN
220	! Only test differences if we are reading in 3 files
221	! Verify time consistency between files
222	IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN
223	CALL ctl_stop( 'Bdy data files must have the same number of time dumps', &
224	& 'Multiple time frequencies not implemented yet' )
225	ENDIF
226	ntimes_bdy = ntimes_bdyt
227	!
228	IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN
229	CALL ctl_stop( 'Bdy data files must have the same time origin', &
230	& 'Multiple time frequencies not implemented yet' )
231	ENDIF
232	zoffset = zoffsett
233	ENDIF
234
235	IF( igrd_start == 1 ) THEN
236	istep(:) = istept(:)
237	ELSE
238	istep(:) = istepu(:)
239	ENDIF
240
241	! Check number of time dumps:
242	IF( ntimes_bdy == 1 .AND. .NOT. ln_bdy_clim ) THEN
243	CALL ctl_stop( 'There is only one time dump in data files', &
244	& 'Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
245	ENDIF
246
247	IF( ln_bdy_clim ) THEN
248	IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN
249	CALL ctl_stop( 'For climatological boundary forcing (ln_bdy_clim=.true.),', &
250	& 'bdy data files must contain 1 or 12 time dumps.' )
251	ELSEIF( ntimes_bdy == 1 ) THEN
252	IF(lwp) WRITE(numout,*)
253	IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
254	ELSEIF( ntimes_bdy == 12 ) THEN
255	IF(lwp) WRITE(numout,*)
256	IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
257	ENDIF
258	ENDIF
259
260	! Find index of first record to read (before first model time).
261	it = 1
262	DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 )
263	it = it + 1
264	END DO
265	nbdy_b = it
266	!
267	WRITE(numout,*) 'Time offset is ',zoffset
268	WRITE(numout,*) 'First record to read is ',nbdy_b
269
270	ENDIF ! endif (nbdy_dta == 1)
271
272
273	! 1.2 Read first record in file if necessary (ie if nbdy_dta == 1)
274	! *****************************************************************
275
276	IF( nbdy_dta == 0) THEN ! boundary data arrays are filled with initial conditions
277	!
278	IF (ln_bdy_tra_frs) THEN
279	igrd = 1 ! T-points data
280	DO ib = 1, nblen(igrd)
281	ii = nbi(ib,igrd)
282	ij = nbj(ib,igrd)
283	DO ik = 1, jpkm1
284	tbdy(ib,ik) = tn(ii, ij, ik)
285	sbdy(ib,ik) = sn(ii, ij, ik)
286	ENDDO
287	END DO
288	ENDIF
289
290	IF(ln_bdy_dyn_frs) THEN
291	igrd = 2 ! U-points data
292	DO ib = 1, nblen(igrd)
293	ii = nbi(ib,igrd)
294	ij = nbj(ib,igrd)
295	DO ik = 1, jpkm1
296	ubdy(ib,ik) = un(ii, ij, ik)
297	ENDDO
298	END DO
299
300	igrd = 3 ! V-points data
301	DO ib = 1, nblen(igrd)
302	ii = nbi(ib,igrd)
303	ij = nbj(ib,igrd)
304	DO ik = 1, jpkm1
305	vbdy(ib,ik) = vn(ii, ij, ik)
306	ENDDO
307	END DO
308	ENDIF
309	!
310	ELSEIF( nbdy_dta == 1 ) THEN ! Set first record in the climatological case:
311	!
312	IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN
313	nbdy_a = 1
314	ELSEIF( ln_bdy_clim .AND. ntimes_bdy == iman ) THEN
315	nbdy_b = 0
316	nbdy_a = imois
317	ELSE
318	nbdy_a = nbdy_b
319	ENDIF
320
321	! Read first record:
322	ipj = 1
323	ipk = jpk
324	igrd = 1
325	ipi = nblendta(igrd)
326
327	IF(ln_bdy_tra_frs) THEN
328	igrd = 1 ! Temperature
329	IF( nblendta(igrd) <= 0 ) THEN
330	idvar = iom_varid( numbdyt, 'votemper' )
331	nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
332	ENDIF
333	WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd)
334	ipi = nblendta(igrd)
335	CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
336
337	DO ib = 1, nblen(igrd)
338	DO ik = 1, jpkm1
339	tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
340	END DO
341	END DO
342	!
343	igrd = 1 ! salinity
344	IF( nblendta(igrd) .le. 0 ) THEN
345	idvar = iom_varid( numbdyt, 'vosaline' )
346	nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
347	ENDIF
348	WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd)
349	ipi = nblendta(igrd)
350	CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
351
352	DO ib = 1, nblen(igrd)
353	DO ik = 1, jpkm1
354	sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
355	END DO
356	END DO
357	ENDIF ! ln_bdy_tra_frs
358
359	IF(ln_bdy_dyn_frs) THEN
360
361	igrd = 2 ! u-velocity
362	IF ( nblendta(igrd) .le. 0 ) THEN
363	idvar = iom_varid( numbdyu,'vozocrtx' )
364	nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar)
365	ENDIF
366	WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd)
367	ipi = nblendta(igrd)
368	CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
369	DO ib = 1, nblen(igrd)
370	DO ik = 1, jpkm1
371	ubdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
372	END DO
373	END DO
374	!
375	igrd = 3 ! v-velocity
376	IF ( nblendta(igrd) .le. 0 ) THEN
377	idvar = iom_varid( numbdyv,'vomecrty' )
378	nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar)
379	ENDIF
380	WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd)
381	ipi = nblendta(igrd)
382	CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
383	DO ib = 1, nblen(igrd)
384	DO ik = 1, jpkm1
385	vbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
386	END DO
387	END DO
388	ENDIF ! ln_bdy_dyn_frs
389
390
391	IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN
392	! First data time is after start of run
393	! Put first value in both time levels
394	nbdy_b = nbdy_a
395	IF(ln_bdy_tra_frs) THEN
396	tbdydta(:,:,1) = tbdydta(:,:,2)
397	sbdydta(:,:,1) = sbdydta(:,:,2)
398	ENDIF
399	IF(ln_bdy_dyn_frs) THEN
400	ubdydta(:,:,1) = ubdydta(:,:,2)
401	vbdydta(:,:,1) = vbdydta(:,:,2)
402	ENDIF
403	END IF
404
405	END IF ! nbdy_dta == 0/1
406
407	! In the case of constant boundary forcing fill bdy arrays once for all
408	IF ((ln_bdy_clim).AND.(ntimes_bdy==1)) THEN
409	IF(ln_bdy_tra_frs) THEN
410	tbdy (:,:) = tbdydta (:,:,2)
411	sbdy (:,:) = sbdydta (:,:,2)
412	ENDIF
413	IF(ln_bdy_dyn_frs) THEN
414	ubdy (:,:) = ubdydta (:,:,2)
415	vbdy (:,:) = vbdydta (:,:,2)
416	ENDIF
417
418	IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )
419	IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu )
420	IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv )
421	END IF
422
423	ENDIF ! End if nit000
424
425
426	! !---------------------!
427	! ! at each time step !
428	! !---------------------!
429
430	IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN
431	!
432	! Read one more record if necessary
433	!**********************************
434
435	IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN ! remember that nbdy_b=0 for kt=nit000
436	nbdy_b = imois
437	nbdy_a = imois + 1
438	nbdy_b = MOD( nbdy_b, iman ) ; IF( nbdy_b == 0 ) nbdy_b = iman
439	nbdy_a = MOD( nbdy_a, iman ) ; IF( nbdy_a == 0 ) nbdy_a = iman
440	lect=.true.
441	ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN
442
443	IF ( nbdy_a < ntimes_bdy ) THEN
444	nbdy_b = nbdy_a
445	nbdy_a = nbdy_a + 1
446	lect =.true.
447	ELSE
448	! We have reached the end of the file
449	! put the last data time into both time levels
450	nbdy_b = nbdy_a
451	IF(ln_bdy_tra_frs) THEN
452	tbdydta(:,:,1) = tbdydta(:,:,2)
453	sbdydta(:,:,1) = sbdydta(:,:,2)
454	ENDIF
455	IF(ln_bdy_dyn_frs) THEN
456	ubdydta(:,:,1) = ubdydta(:,:,2)
457	vbdydta(:,:,1) = vbdydta(:,:,2)
458	ENDIF
459	END IF ! nbdy_a < ntimes_bdy
460
461	END IF
462
463	IF( lect ) THEN
464	! Swap arrays
465	IF(ln_bdy_tra_frs) THEN
466	tbdydta(:,:,1) = tbdydta(:,:,2)
467	sbdydta(:,:,1) = sbdydta(:,:,2)
468	ENDIF
469	IF(ln_bdy_dyn_frs) THEN
470	ubdydta(:,:,1) = ubdydta(:,:,2)
471	vbdydta(:,:,1) = vbdydta(:,:,2)
472	ENDIF
473
474	! read another set
475	ipj = 1
476	ipk = jpk
477
478	IF(ln_bdy_tra_frs) THEN
479	!
480	igrd = 1 ! temperature
481	ipi = nblendta(igrd)
482	CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
483	DO ib = 1, nblen(igrd)
484	DO ik = 1, jpkm1
485	tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
486	END DO
487	END DO
488	!
489	igrd = 1 ! salinity
490	ipi = nblendta(igrd)
491	CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
492	DO ib = 1, nblen(igrd)
493	DO ik = 1, jpkm1
494	sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
495	END DO
496	END DO
497	ENDIF ! ln_bdy_tra_frs
498
499	IF(ln_bdy_dyn_frs) THEN
500	!
501	igrd = 2 ! u-velocity
502	ipi = nblendta(igrd)
503	CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
504	DO ib = 1, nblen(igrd)
505	DO ik = 1, jpkm1
506	ubdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
507	END DO
508	END DO
509	!
510	igrd = 3 ! v-velocity
511	ipi = nblendta(igrd)
512	CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
513	DO ib = 1, nblen(igrd)
514	DO ik = 1, jpkm1
515	vbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
516	END DO
517	END DO
518	ENDIF ! ln_bdy_dyn_frs
519
520	!
521	IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b
522	IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy_a ',nbdy_a
523	IF (.NOT.ln_bdy_clim) THEN
524	IF(lwp) WRITE(numout,*) 'first record time (s): ', istep(nbdy_b)
525	IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer
526	IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a)
527	ENDIF
528	!
529	ENDIF ! end lect=.true.
530
531
532	! Interpolate linearly
533	! ********************
534	!
535	IF( ln_bdy_clim ) THEN ; zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b) ) + 0.5 - i15
536	ELSE ; zxy = FLOAT( istep(nbdy_b) - itimer ) / FLOAT( istep(nbdy_b) - istep(nbdy_a) )
537	END IF
538
539	IF(ln_bdy_tra_frs) THEN
540	igrd = 1 ! temperature & salinity
541	DO ib = 1, nblen(igrd)
542	DO ik = 1, jpkm1
543	tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1)
544	sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1)
545	END DO
546	END DO
547	ENDIF
548
549	IF(ln_bdy_dyn_frs) THEN
550	igrd = 2 ! u-velocity
551	DO ib = 1, nblen(igrd)
552	DO ik = 1, jpkm1
553	ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1)
554	END DO
555	END DO
556	!
557	igrd = 3 ! v-velocity
558	DO ib = 1, nblen(igrd)
559	DO ik = 1, jpkm1
560	vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1)
561	END DO
562	END DO
563	ENDIF
564
565	END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy>1))
566
567
568	! !---------------------!
569	! ! last call !
570	! !---------------------!
571	IF( kt == nitend ) THEN
572	IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt ) ! Closing of the 3 files
573	IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu )
574	IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv )
575	ENDIF
576	!
577	ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs
578
579	END SUBROUTINE bdy_dta
580
581
582	SUBROUTINE bdy_dta_bt( kt, jit )
583	!!---------------------------------------------------------------------------
584	!! * SUBROUTINE bdy_dta_bt *
585	!!
586	!! ** Purpose : Read unstructured boundary data for Flather condition
587	!!
588	!! ** Method : At the first timestep, read in boundary data for two
589	!! times from the file and time-interpolate. At other
590	!! timesteps, check to see if we need another time from
591	!! the file. If so read it in. Time interpolate.
592	!!---------------------------------------------------------------------------
593	!!gm DOCTOR names : argument integer : start with "k"
594	INTEGER, INTENT( in ) :: kt ! ocean time-step index
595	INTEGER, INTENT( in ) :: jit ! barotropic time step index
596	! ! (for timesplitting option, otherwise zero)
597	!!
598	LOGICAL :: lect ! flag for reading
599	INTEGER :: it, ib, igrd ! dummy loop indices
600	INTEGER :: idvar ! netcdf var ID
601	INTEGER :: iman, i15, imois ! Time variables for monthly clim forcing
602	INTEGER :: ntimes_bdyt, ntimes_bdyu, ntimes_bdyv
603	INTEGER :: itimer, totime
604	INTEGER :: ipi, ipj, ipk, inum ! temporary integers (NetCDF read)
605	INTEGER :: iyear0, imonth0, iday0
606	INTEGER :: ihours0, iminutes0, isec0
607	INTEGER :: iyear, imonth, iday, isecs
608	INTEGER, DIMENSION(jpbtime) :: istept, istepu, istepv ! time arrays from data files
609	REAL(wp) :: zxy, zoffsett
610	REAL(wp) :: zoffsetu, zoffsetv
611	REAL(wp) :: dayjul0, zdayjulini
612	REAL(wp) :: zinterval_s, zinterval_e ! First and last interval in time axis
613	REAL(wp), DIMENSION(jpbtime) :: zstepr ! REAL time array from data files
614	REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields
615	CHARACTER(LEN=80), DIMENSION(3) :: clfile
616	!!---------------------------------------------------------------------------
617
618	!!gm add here the same style as in bdy_dta
619	!!gm clearly bdy_dta_bt and bdy_dta can be combined...
620	!!gm too many things duplicated in the read of data... simplification can be done
621
622	! -------------------- !
623	! Initialization !
624	! -------------------- !
625
626	lect = .false. ! If true, read a time record
627
628	! Some time variables for monthly climatological forcing:
629	! *******************************************************
630	!!gm here use directely daymod variables
631
632	iman = INT( raamo ) ! Number of months in a year
633
634	i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
635	! i15=0 if the current day is in the first half of the month, else i15=1
636
637	imois = nmonth + i15 - 1 ! imois is the first month record
638	IF( imois == 0 ) imois = iman
639
640	! Time variable for non-climatological forcing:
641	! *********************************************
642
643	itimer = ((kt-1)-nit000+1)*rdt ! current time in seconds for interpolation
644	itimer = itimer + jit*rdt/REAL(nn_baro,wp) ! in non-climatological case
645
646	IF ( ln_bdy_tides ) THEN
647
648	! -------------------------------------!
649	! Update BDY fields with tidal forcing !
650	! -------------------------------------!
651
652	CALL tide_update( kt, jit )
653
654	ENDIF
655
656	IF ( ln_bdy_dyn_fla ) THEN
657
658	! -------------------------------------!
659	! Update BDY fields with model data !
660	! -------------------------------------!
661
662	! !-------------------!
663	IF( kt == nit000 ) THEN ! First call only !
664	! !-------------------!
665	istep_bt(:) = 0
666	nbdy_b_bt = 0
667	nbdy_a_bt = 0
668
669	! Get time information from bdy data file
670	! ***************************************
671
672	IF(lwp) WRITE(numout,*)
673	IF(lwp) WRITE(numout,*) 'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.'
674	IF(lwp) WRITE(numout,*) '~~~~~~~'
675
676	IF( nbdy_dta == 0 ) THEN
677	IF(lwp) WRITE(numout,*) 'Bdy data are taken from initial conditions'
678
679	ELSEIF (nbdy_dta == 1) THEN
680	IF(lwp) WRITE(numout,*) 'Bdy data are read in netcdf files'
681
682	totime = (nitend-nit000+1)*rdt ! Total time of the run to verify that all the
683	! necessary time dumps in file are included
684
685	clfile(1) = filbdy_data_bt_T
686	clfile(2) = filbdy_data_bt_U
687	clfile(3) = filbdy_data_bt_V
688
689	DO igrd = 1,3
690
691	CALL iom_open( clfile(igrd), inum )
692	CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, offset=zoffset, &
693	zyear=nyear, zmonth=nmonth, zday=nday, zsec_day=rsec_day )
694
695	SELECT CASE( igrd )
696	CASE (1)
697	numbdyt = inum
698	CASE (2)
699	numbdyu = inum
700	CASE (3)
701	numbdyv = inum
702	END SELECT
703
704	!! TO BE DONE... Check consistency between calendar from file
705	!! (available optionally from iom_gettime) and calendar in model
706	!! when calendar in model available outside of IOIPSL.
707
708	! Check that there are not too many times in the file.
709	IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( &
710	'Number of time dumps in bdy file exceed jpbtime parameter', &
711	'Check file:' // TRIM(clfile(igrd)) )
712
713	! Check that time array increases (or interp will fail):
714	DO it = 2, ntimes_bdy
715	IF ( zstepr(it-1) >= zstepr(it) ) THEN
716	CALL ctl_stop('Time array in unstructured boundary data file', &
717	'does not continuously increase.', &
718	'Check file:' // TRIM(clfile(igrd)) )
719	EXIT
720	END IF
721	END DO
722
723	IF ( .NOT. ln_bdy_clim ) THEN
724	! Check that times in file span model run time:
725
726	! Note: the fields may be time means, so we allow nit000 to be before
727	! first time in the file, provided that it falls inside the meaning
728	! period of the first field. Until we can get the meaning period
729	! from the file, use the interval between fields as a proxy.
730	! If nit000 is before the first time, use the value at first time
731	! instead of extrapolating. This is done by putting time 1 into
732	! both time levels.
733	! The same applies to the last time level: see setting of lect below.
734
735	IF ( ntimes_bdy == 1 ) CALL ctl_stop( &
736	'There is only one time dump in data files', &
737	'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
738
739	zinterval_s = zstepr(2) - zstepr(1)
740	zinterval_e = zstepr(ntimes_bdy) - zstepr(ntimes_bdy-1)
741
742	IF ( zstepr(1) - zinterval_s / 2.0 > 0 ) THEN
743	IF(lwp) WRITE(numout,*) 'First bdy time relative to nit000:', zstepr(1)
744	IF(lwp) WRITE(numout,*) 'Interval between first two times: ', zinterval_s
745	CALL ctl_stop( 'First data time is after start of run', &
746	'by more than half a meaning period', &
747	'Check file: ' // TRIM(clfile(igrd)) )
748	END IF
749
750	IF ( zstepr(ntimes_bdy) + zinterval_e / 2.0 < totime ) THEN
751	IF(lwp) WRITE(numout,*) 'Last bdy time relative to nit000:', zstepr(ntimes_bdy)
752	IF(lwp) WRITE(numout,*) 'Interval between last two times: ', zinterval_e
753	CALL ctl_stop( 'Last data time is before end of run', &
754	'by more than half a meaning period', &
755	'Check file: ' // TRIM(clfile(igrd)) )
756	END IF
757
758	END IF ! .NOT. ln_bdy_clim
759
760	IF ( igrd .EQ. 1) THEN
761	ntimes_bdyt = ntimes_bdy_bt
762	zoffsett = zoffset
763	istept(:) = INT( zstepr(:) + zoffset )
764	ELSE IF (igrd .EQ. 2) THEN
765	ntimes_bdyu = ntimes_bdy_bt
766	zoffsetu = zoffset
767	istepu(:) = INT( zstepr(:) + zoffset )
768	ELSE IF (igrd .EQ. 3) THEN
769	ntimes_bdyv = ntimes_bdy_bt
770	zoffsetv = zoffset
771	istepv(:) = INT( zstepr(:) + zoffset )
772	ENDIF
773
774	ENDDO
775
776	! Verify time consistency between files
777
778	IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN
779	CALL ctl_stop( &
780	'Time axis lengths differ between bdy data files', &
781	'Multiple time frequencies not implemented yet' )
782	ELSE
783	ntimes_bdy_bt = ntimes_bdyt
784	ENDIF
785
786	IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN
787	CALL ctl_stop( &
788	'Bdy data files must have the same time origin', &
789	'Multiple time frequencies not implemented yet' )
790	ENDIF
791	zoffset = zoffsett
792
793	!! Check that times are the same in the three files... HERE.
794	istep_bt(:) = istept(:)
795
796	! Check number of time dumps:
797	IF (ln_bdy_clim) THEN
798	SELECT CASE ( ntimes_bdy_bt )
799	CASE( 1 )
800	IF(lwp) WRITE(numout,*)
801	IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
802	IF(lwp) WRITE(numout,*)
803	CASE( 12 )
804	IF(lwp) WRITE(numout,*)
805	IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
806	IF(lwp) WRITE(numout,*)
807	CASE DEFAULT
808	CALL ctl_stop( &
809	'For climatological boundary forcing (ln_bdy_clim=.true.),',&
810	'bdy data files must contain 1 or 12 time dumps.' )
811	END SELECT
812	ENDIF
813
814	! Find index of first record to read (before first model time).
815
816	it=1
817	DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1)))
818	it=it+1
819	END DO
820	nbdy_b_bt = it
821
822	WRITE(numout,*) 'Time offset is ',zoffset
823	WRITE(numout,*) 'First record to read is ',nbdy_b_bt
824
825	ENDIF ! endif (nbdy_dta == 1)
826
827	! 1.2 Read first record in file if necessary (ie if nbdy_dta == 1)
828	! *****************************************************************
829
830	IF ( nbdy_dta == 0) THEN
831	! boundary data arrays are filled with initial conditions
832	igrd = 2 ! U-points data
833	DO ib = 1, nblen(igrd)
834	ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1)
835	END DO
836
837	igrd = 3 ! V-points data
838	DO ib = 1, nblen(igrd)
839	vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1)
840	END DO
841
842	igrd = 1 ! T-points data
843	DO ib = 1, nblen(igrd)
844	sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd))
845	END DO
846
847	ELSEIF (nbdy_dta == 1) THEN
848
849	! Set first record in the climatological case:
850	IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN
851	nbdy_a_bt = 1
852	ELSEIF ((ln_bdy_clim).AND.(ntimes_bdy_bt==iman)) THEN
853	nbdy_b_bt = 0
854	nbdy_a_bt = imois
855	ELSE
856	nbdy_a_bt = nbdy_b_bt
857	END IF
858
859	! Open Netcdf files:
860
861	CALL iom_open ( filbdy_data_bt_T, numbdyt_bt )
862	CALL iom_open ( filbdy_data_bt_U, numbdyu_bt )
863	CALL iom_open ( filbdy_data_bt_V, numbdyv_bt )
864
865	! Read first record:
866	ipj=1
867	igrd=1
868	ipi=nblendta(igrd)
869
870	! ssh
871	igrd=1
872	IF ( nblendta(igrd) .le. 0 ) THEN
873	idvar = iom_varid( numbdyt_bt,'sossheig' )
874	nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar)
875	ENDIF
876	WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd)
877	ipi=nblendta(igrd)
878
879	CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt )
880
881	DO ib=1, nblen(igrd)
882	sshbdydta(ib,2) = zdta(nbmap(ib,igrd),1)
883	END DO
884
885	! u-velocity
886	igrd=2
887	IF ( nblendta(igrd) .le. 0 ) THEN
888	idvar = iom_varid( numbdyu_bt,'vobtcrtx' )
889	nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar)
890	ENDIF
891	WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd)
892	ipi=nblendta(igrd)
893
894	CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt )
895
896	DO ib=1, nblen(igrd)
897	ubtbdydta(ib,2) = zdta(nbmap(ib,igrd),1)
898	END DO
899
900	! v-velocity
901	igrd=3
902	IF ( nblendta(igrd) .le. 0 ) THEN
903	idvar = iom_varid( numbdyv_bt,'vobtcrty' )
904	nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar)
905	ENDIF
906	WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd)
907	ipi=nblendta(igrd)
908
909	CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt )
910
911	DO ib=1, nblen(igrd)
912	vbtbdydta(ib,2) = zdta(nbmap(ib,igrd),1)
913	END DO
914
915	END IF
916
917	! In the case of constant boundary forcing fill bdy arrays once for all
918	IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN
919
920	ubtbdy (:) = ubtbdydta (:,2)
921	vbtbdy (:) = vbtbdydta (:,2)
922	sshbdy (:) = sshbdydta (:,2)
923
924	CALL iom_close( numbdyt_bt )
925	CALL iom_close( numbdyu_bt )
926	CALL iom_close( numbdyv_bt )
927
928	END IF
929
930	ENDIF ! End if nit000
931
932	! -------------------- !
933	! 2. At each time step !
934	! -------------------- !
935
936	IF ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) THEN
937
938	! 2.1 Read one more record if necessary
939	!**************************************
940
941	IF ( (ln_bdy_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000
942	nbdy_b_bt = imois
943	nbdy_a_bt = imois+1
944	nbdy_b_bt = MOD( nbdy_b_bt, iman )
945	IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman
946	nbdy_a_bt = MOD( nbdy_a_bt, iman )
947	IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman
948	lect=.true.
949
950	ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN
951	nbdy_b_bt=nbdy_a_bt
952	nbdy_a_bt=nbdy_a_bt+1
953	lect=.true.
954	END IF
955
956	IF (lect) THEN
957
958	! Swap arrays
959	sshbdydta(:,1) = sshbdydta(:,2)
960	ubtbdydta(:,1) = ubtbdydta(:,2)
961	vbtbdydta(:,1) = vbtbdydta(:,2)
962
963	! read another set
964
965	ipj=1
966	ipk=jpk
967	igrd=1
968	ipi=nblendta(igrd)
969
970
971	! ssh
972	igrd=1
973	ipi=nblendta(igrd)
974
975	CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt )
976
977	DO ib=1, nblen(igrd)
978	sshbdydta(ib,2) = zdta(nbmap(ib,igrd),1)
979	END DO
980
981	! u-velocity
982	igrd=2
983	ipi=nblendta(igrd)
984
985	CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt )
986
987	DO ib=1, nblen(igrd)
988	ubtbdydta(ib,2) = zdta(nbmap(ib,igrd),1)
989	END DO
990
991	! v-velocity
992	igrd=3
993	ipi=nblendta(igrd)
994
995	CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt )
996
997	DO ib=1, nblen(igrd)
998	vbtbdydta(ib,2) = zdta(nbmap(ib,igrd),1)
999	END DO
1000
1001
1002	IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b_bt ',nbdy_b_bt
1003	IF(lwp) WRITE(numout,*) '~~~~~~~~ last record file used nbdy_a_bt ',nbdy_a_bt
1004	IF (.NOT.ln_bdy_clim) THEN
1005	IF(lwp) WRITE(numout,*) 'first record time (s): ', istep_bt(nbdy_b_bt)
1006	IF(lwp) WRITE(numout,*) 'model time (s) : ', itimer
1007	IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt)
1008	ENDIF
1009	END IF ! end lect=.true.
1010
1011
1012	! 2.2 Interpolate linearly:
1013	! ***************************
1014
1015	IF (ln_bdy_clim) THEN
1016	zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b_bt) ) + 0.5 - i15
1017	ELSE
1018	zxy = FLOAT(istep_bt(nbdy_b_bt)-itimer) / FLOAT(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt))
1019	END IF
1020
1021	igrd=1
1022	DO ib=1, nblen(igrd)
1023	sshbdy(ib) = zxy * sshbdydta(ib,2) + &
1024	(1.-zxy) * sshbdydta(ib,1)
1025	END DO
1026
1027	igrd=2
1028	DO ib=1, nblen(igrd)
1029	ubtbdy(ib) = zxy * ubtbdydta(ib,2) + &
1030	(1.-zxy) * ubtbdydta(ib,1)
1031	END DO
1032
1033	igrd=3
1034	DO ib=1, nblen(igrd)
1035	vbtbdy(ib) = zxy * vbtbdydta(ib,2) + &
1036	(1.-zxy) * vbtbdydta(ib,1)
1037	END DO
1038
1039
1040	END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy_bt>1))
1041
1042	! ------------------- !
1043	! Last call kt=nitend !
1044	! ------------------- !
1045
1046	! Closing of the 3 files
1047	IF( kt == nitend ) THEN
1048	CALL iom_close( numbdyt_bt )
1049	CALL iom_close( numbdyu_bt )
1050	CALL iom_close( numbdyv_bt )
1051	ENDIF
1052
1053	ENDIF ! ln_bdy_dyn_frs
1054
1055	END SUBROUTINE bdy_dta_bt
1056
1057
1058	#else
1059	!!----------------------------------------------------------------------
1060	!! Dummy module NO Unstruct Open Boundary Conditions
1061	!!----------------------------------------------------------------------
1062	CONTAINS
1063	SUBROUTINE bdy_dta( kt ) ! Empty routine
1064	WRITE(,) 'bdy_dta: You should not have seen this print! error?', kt
1065	END SUBROUTINE bdy_dta
1066	SUBROUTINE bdy_dta_bt( kt, kit ) ! Empty routine
1067	WRITE(,) 'bdy_dta: You should not have seen this print! error?', kt, kit
1068	END SUBROUTINE bdy_dta_bt
1069	#endif
1070
1071	!!==============================================================================
1072	END MODULE bdydta

Download in other formats:

Original Format

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

Ticket #377: bdydta.F90

Download in other formats: