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.

bdydta.F90 in trunk/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/BDY/bdydta.F90 @ 1110

Last change on this file since 1110 was 1058, checked in by rblod, 16 years ago
Correct preprocessing syntax, see ticket #160
Property svn:executable set to ``*
File size: 39.3 KB

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

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

Download in other formats: