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

Last change on this file since 35 was 35, checked in by opalod, 20 years ago
CT : UPDATE002 : Concerns Open Boundaries
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.6 KB

Line
1	MODULE obcdta
2	!!==============================================================================
3	!! * MODULE obcdta *
4	!! Open boundary data : read the data for the open boundaries.
5	!!==============================================================================
6	#if defined key_obc
7	!!------------------------------------------------------------------------------
8	!! 'key_obc' : Open Boundary Conditions
9	!!------------------------------------------------------------------------------
10	!! obc_dta : read u, v, t, s data along each open boundary
11	!! obc_dta_psi : read psi data along each open boundary (rigid lid only)
12	!!------------------------------------------------------------------------------
13	!! * Modules used
14	USE oce ! ocean dynamics and tracers
15	USE dom_oce ! ocean space and time domain
16	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
17	USE phycst ! physical constants
18	USE obc_oce ! ocean open boundary conditions
19	USE daymod ! calendar
20	USE in_out_manager ! I/O logical units
21	USE lib_mpp ! distribued memory computing
22	USE dynspg_rl !
23
24
25	# if ! defined key_dynspg_fsc
26	USE obccli
27	# endif
28
29	IMPLICIT NONE
30	PRIVATE
31
32	!! * Accessibility
33	PUBLIC obc_dta ! routines called by step.F90
34
35	!! * Substitutions
36	# include "obc_vectopt_loop_substitute.h90"
37	!!---------------------------------------------------------------------------------
38	!! OPA 9.0 , LODYC-IPSL (2003)
39	!!---------------------------------------------------------------------------------
40
41	CONTAINS
42
43	SUBROUTINE obc_dta( kt )
44	!!---------------------------------------------------------------------------
45	!! * SUBROUTINE obc_dta *
46	!!
47	!! ** Purpose : Find the climatological boundary arrays for the specified date,
48	!! Originally this routine interpolated between monthly fields
49	!! of a climatology.
50	!! When using hydrological section data, we have only on snapshot
51	!! and do not need to interpolate.
52	!!
53	!! ** Method : Determine the current month from kdat, and interpolate for
54	!! the current day.
55	!!
56	!! History :
57	!! ! 98-05 (J.M. Molines) Original code
58	!! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
59	!!---------------------------------------------------------------------------
60	!! * Arguments
61	INTEGER, INTENT( in ) :: kt ! ocean time-step index
62
63	!! * Local declarations
64	INTEGER :: ji, jj, jk, ii, ij ! dummy loop indices
65	INTEGER :: inum = 11 ! temporary logical unit
66	INTEGER :: imois, iman, imoisp1
67	INTEGER :: i15
68	INTEGER :: irecl, irec, ios
69	INTEGER, PARAMETER :: &
70	jpmois = 1, &
71	jpf = 3
72	REAL(wp) :: zxy
73	CHARACTER (len=4 ) :: clversion
74	CHARACTER (len=80) :: clcom
75	!!---------------------------------------------------------------------
76
77
78	IF( lk_dynspg_rl ) CALL obc_dta_psi( kt ) ! update bsf data at open boundaries
79
80
81	! 0. Initialization of date
82	! imois is the index (1 to 12) of the first month to be used in the
83	! interpolation. ndastp is the date (integer format yymmdd) calculated
84	! in routine day.F starting from ndate0 given in namelist.
85	! -----------------------------------------------------------------------
86
87	iman = jpmois
88	i15 = nday/16
89	imois = nmonth + i15 - 1
90	IF( imois == 0 ) imois = iman
91	imoisp1 = MOD( imois + 1, iman )
92	IF( imoisp1 == 0 ) imoisp1 = iman
93	! ... zxy is the weight of the after field
94	zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
95	IF( zxy > 1.01 .OR. zxy < 0. ) THEN
96	IF(lwp) WRITE(numout,*)' '
97	IF(lwp) WRITE(numout,*)'obc_dta: Pbm with the the weight of the after field zxy '
98	IF(lwp) WRITE(numout,*)'~~~~~~~~~~~'
99	nstop = nstop + 1
100	END IF
101
102	! 1. First call:
103	! ----------------
104
105	IF( kt == nit000 ) THEN
106	IF(lwp) WRITE(numout,*)' '
107	IF(lwp) WRITE(numout,*)'obcdta: initial step in obc_dta'
108	IF(lwp) WRITE(numout,*)'~~~~~~ months ',imois,' and', imoisp1,' read'
109
110	! 1.1 Tangential velocities set to zero
111	! --------------------------------------
112	IF( lpeastobc ) THEN
113	DO jk=1, jpkm1
114	vfoe(:,jk)=0.00
115	END DO
116	END IF
117
118	IF( lpwestobc ) THEN
119	DO jk=1, jpkm1
120	vfow(:,jk)=0.00
121	END DO
122	END IF
123
124	IF( lpsouthobc ) THEN
125	DO jk=1, jpkm1
126	ufos(:,jk)=0.00
127	END DO
128	END IF
129
130	IF( lpnorthobc ) THEN
131	DO jk=1, jpkm1
132	ufon(:,jk)=0.00
133	END DO
134	END IF
135
136	! 1.2 Set the Normal velocity and tracers data for the EAST OBC
137	! -------------------------------------------------------------
138
139	IF( lpeastobc ) THEN
140
141	! initialisation to zero
142	sedta(:,:,1) = 0.e0
143	tedta(:,:,1) = 0.e0
144	uedta(:,:,1) = 0.e0
145	! ! ================== !
146	IF( nobc_dta == 0 ) THEN ! initial state used
147	! ! ================== !
148	DO ji = fs_nie0, fs_nie1 ! Vector opt.
149	DO jk = 1, jpkm1
150	DO jj = 1, jpj
151	ij = jj -1 + njmpp
152	sedta(ij,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk)
153	tedta(ij,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk)
154	uedta(ij,jk,1) = un(ji,jj,jk)*umask(ji,jj,jk)
155	END DO
156	END DO
157	END DO
158
159	DO jk = 1, jpkm1
160	DO jj = 1, jpj
161	ij = jj -1 + njmpp
162	sfoe(jj,jk) = sedta(ij,jk,1)*temsk(jj,jk)
163	tfoe(jj,jk) = tedta(ij,jk,1)*temsk(jj,jk)
164	ufoe(jj,jk) = uedta(ij,jk,1)*uemsk(jj,jk)
165	END DO
166	END DO
167	! ! =================== !
168	ELSE ! read in obceast.dta
169	! ! =================== !
170	OPEN(UNIT = inum, &
171	IOSTAT = ios, &
172	FILE ='obceast.dta', &
173	FORM ='UNFORMATTED', &
174	ACCESS ='DIRECT', &
175	RECL = 4096 )
176	IF( ios > 0 ) THEN
177	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obceast.dta file '
178	IF(lwp) WRITE(numout,*) '~~~~~~~'
179	nstop = nstop + 1
180	END IF
181	READ(inum,REC=1) clversion,clcom,irecl
182	CLOSE(inum)
183	IF(lwp) WRITE(numout,*)' '
184	IF(lwp) WRITE(numout,*)' opening obceast.dta with irecl=',irecl
185	OPEN(UNIT = inum, &
186	IOSTAT = ios, &
187	FILE ='obceast.dta', &
188	FORM ='UNFORMATTED', &
189	ACCESS ='DIRECT', &
190	RECL = irecl )
191	IF( ios > 0 ) THEN
192	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obceast.dta file '
193	IF(lwp) WRITE(numout,*) '~~~~~~~'
194	nstop = nstop + 1
195	END IF
196
197	! ... Read datafile and set temperature, salinity and normal velocity
198	! ... initialise the sedta, tedta arrays
199	! ... index 1 refer to before, 2 to after
200
201	DO jk = 1, jpkm1
202	irec = 2 + (jk -1)* jpf
203	READ(inum,REC=irec )((sedta(jj,jk,1),ji=1,1),jj=jpjed, jpjef)
204	READ(inum,REC=irec+1)((tedta(jj,jk,1),ji=1,1),jj=jpjed, jpjef)
205	READ(inum,REC=irec+2)((uedta(jj,jk,1),ji=1,1),jj=jpjed, jpjef)
206	! ... set climatological temperature and salinity at the boundary
207	! ... do not interpolate between months : impose values of climatology
208	DO jj = 1, jpj
209	ij = jj -1 + njmpp
210	sfoe(jj,jk) = sedta(ij,jk,1)*temsk(jj,jk)
211	tfoe(jj,jk) = tedta(ij,jk,1)*temsk(jj,jk)
212	END DO
213	END DO
214	CLOSE(inum)
215
216	# if ! defined key_dynspg_fsc
217	! ... Rigid lid case: make sure uedta is baroclinic velocity
218	! ... In rigid lid case uedta needs to be the baroclinic component.
219
220	CALL obc_cli( uedta, uclie, fs_nie0, fs_nie1, 0, jpj, njmpp )
221
222	# endif
223	! ... Set normal velocity (on nie0, nie1 <=> jpieob)
224	DO jk = 1, jpkm1
225	DO jj = 1, jpj
226	ij = jj -1 + njmpp
227	ufoe(jj,jk) = uedta(ij,jk,1)*uemsk(jj,jk)
228	END DO
229	END DO
230	ENDIF
231	ENDIF
232
233	! 1.3 Set the Normal velocity and tracers data for the WEST OBC
234	! -------------------------------------------------------------
235
236	IF( lpwestobc) THEN
237
238	! initialisation to zero
239	swdta(:,:,1) = 0.e0
240	twdta(:,:,1) = 0.e0
241	uwdta(:,:,1) = 0.e0
242	! ! ================== !
243	IF( nobc_dta == 0 ) THEN ! initial state used
244	! ! ================== !
245	DO ji = fs_niw0, fs_niw1 ! Vector opt.
246	DO jk = 1, jpkm1
247	DO jj = 1, jpj
248	ij = jj -1 + njmpp
249	swdta(ij,jk,1) = sn(ji,jj,jk)*tmask(ji,jj,jk)
250	twdta(ij,jk,1) = tn(ji,jj,jk)*tmask(ji,jj,jk)
251	uwdta(ij,jk,1) = un(ji,jj,jk)*umask(ji,jj,jk)
252	END DO
253	END DO
254	END DO
255
256	DO jk = 1, jpkm1
257	DO jj = 1, jpj
258	ij = jj -1 + njmpp
259	sfow(jj,jk) = swdta(ij,jk,1)*twmsk(jj,jk)
260	tfow(jj,jk) = twdta(ij,jk,1)*twmsk(jj,jk)
261	ufow(jj,jk) = uwdta(ij,jk,1)*uwmsk(jj,jk)
262	END DO
263	END DO
264	! ! =================== !
265	ELSE ! read in obceast.dta
266	! ! =================== !
267	OPEN(UNIT = inum, &
268	IOSTAT = ios, &
269	FILE ='obcwest.dta', &
270	FORM ='UNFORMATTED', &
271	ACCESS ='DIRECT', &
272	RECL = 4096 )
273	IF( ios > 0 ) THEN
274	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obcwest.dta file '
275	IF(lwp) WRITE(numout,*) '~~~~~~~'
276	nstop = nstop + 1
277	END IF
278	READ(inum,REC=1) clversion, clcom,irecl
279	CLOSE(inum)
280	IF(lwp) WRITE(numout,*)' '
281	IF(lwp) WRITE(numout,*)' opening obcwest.dta with irecl=',irecl
282	OPEN(UNIT = inum, &
283	IOSTAT = ios, &
284	FILE ='obcwest.dta', &
285	FORM ='UNFORMATTED', &
286	ACCESS ='DIRECT', &
287	RECL = irecl )
288	IF( ios > 0 ) THEN
289	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obcwest.dta file '
290	IF(lwp) WRITE(numout,*) '~~~~~~~'
291	nstop = nstop + 1
292	END IF
293
294	! ... Read datafile and set temperature, salinity and normal velocity
295	! ... initialise the swdta, twdta arrays
296	! ... index 1 refer to before, 2 to after
297	DO jk = 1, jpkm1
298	irec = 2 + (jk -1)* jpf
299	READ(inum,REC=irec )((swdta(jj,jk,1),ji=1,1),jj=jpjwd, jpjwf)
300	READ(inum,REC=irec+1)((twdta(jj,jk,1),ji=1,1),jj=jpjwd, jpjwf)
301	READ(inum,REC=irec+2)((uwdta(jj,jk,1),ji=1,1),jj=jpjwd, jpjwf)
302	DO jj = 1, jpj
303	ij = jj -1 + njmpp
304	sfow(jj,jk) = swdta(ij,jk,1)*twmsk(jj,jk)
305	tfow(jj,jk) = twdta(ij,jk,1)*twmsk(jj,jk)
306	END DO
307	END DO
308	CLOSE(inum)
309
310	#if ! defined key_dynspg_fsc
311	! ... Rigid lid case: make sure uwdta is baroclinic velocity
312	! ... In rigid lid case uwdta needs to be the baroclinic component.
313
314	CALL obc_cli( uwdta, ucliw, fs_niw0, fs_niw1, 0, jpj, njmpp )
315
316	# endif
317	! ... Set normal velocity (on niw0, niw1 <=> jpiwob)
318	DO jk = 1, jpkm1
319	DO jj = 1, jpj
320	ij = jj -1 + njmpp
321	ufow(jj,jk) = uwdta(ij,jk,1)*uwmsk(jj,jk)
322	END DO
323	END DO
324	ENDIF
325	ENDIF
326
327	! 1.4 Set the Normal velocity and tracers data for the NORTH OBC
328	! ---------------------------------------------------------------
329
330	IF( lpnorthobc) THEN
331
332	! initialisation to zero
333	sndta(:,:,1) = 0.e0
334	tndta(:,:,1) = 0.e0
335	vndta(:,:,1) = 0.e0
336
337	OPEN(UNIT = inum, &
338	IOSTAT = ios, &
339	FILE ='obcnorth.dta', &
340	FORM ='UNFORMATTED', &
341	ACCESS ='DIRECT', &
342	RECL = 4096 )
343	IF( ios > 0 ) THEN
344	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obcnorth.dta file '
345	IF(lwp) WRITE(numout,*) '~~~~~~~'
346	nstop = nstop + 1
347	END IF
348	READ(inum,REC=1) clversion,clcom,irecl
349	CLOSE(inum)
350	IF(lwp) WRITE(numout,*)' '
351	IF(lwp) WRITE(numout,*)' opening obcnorth.dta with irecl=',irecl
352	OPEN(UNIT = inum, &
353	IOSTAT = ios, &
354	FILE ='obcnorth.dta', &
355	FORM ='UNFORMATTED', &
356	ACCESS ='DIRECT', &
357	RECL = irecl )
358	IF( ios > 0 ) THEN
359	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obcnorth.dta file '
360	IF(lwp) WRITE(numout,*) '~~~~~~~'
361	nstop = nstop + 1
362	END IF
363
364	! ... Read datafile and set temperature and salinity
365	! ... initialise the sndta, tndta arrays
366	! ... index 1 refer to before, 2 to after
367	DO jk = 1, jpkm1
368	irec = 2 + (jk -1)* jpf
369	READ(inum,REC=irec )((sndta(jj,jk,1),ji=1,1),jj=jpind, jpinf)
370	READ(inum,REC=irec+1)((tndta(jj,jk,1),ji=1,1),jj=jpind, jpinf)
371	READ(inum,REC=irec+2)((vndta(jj,jk,1),ji=1,1),jj=jpind, jpinf)
372	! ... do not interpolate: impose values of climatology
373	DO ji = 1, jpi
374	ii = ji -1 + nimpp
375	sfon(ji,jk) = sndta(ii,jk,1)*tnmsk(ji,jk)
376	tfon(ji,jk) = tndta(ii,jk,1)*tnmsk(ji,jk)
377	END DO
378	END DO
379	CLOSE(inum)
380
381	# if ! defined key_dynspg_fsc
382	! ... Rigid lid case: make sure vndta is baroclinic velocity
383	! In rigid lid case vndta needs to be the baroclinic component.
384	! substract the barotropic velocity component (zvnbtpe)
385	! from the total one (vndta).
386
387	CALL obc_cli( vndta, vclin, fs_njn0, fs_njn1, 1, jpi, nimpp )
388	# endif
389	! ... Set normal velocity (on njn0, njn1 <=> jpjnob)
390	DO jk = 1, jpkm1
391	DO ji = 1, jpi
392	ii = ji -1 + nimpp
393	vfon(ji,jk) = vndta(ii,jk,1)*vnmsk(ji,jk)
394	END DO
395	END DO
396
397	END IF
398
399	! 1.5 Set the Normal velocity and tracers data for the SOUTH OBC
400	! ---------------------------------------------------------------
401
402	IF( lpsouthobc) THEN
403
404	! initialisation to zero
405	ssdta(:,:,1) = 0.e0
406	tsdta(:,:,1) = 0.e0
407	vsdta(:,:,1) = 0.e0
408
409	OPEN(UNIT = inum, &
410	IOSTAT = ios, &
411	FILE ='obcsouth.dta', &
412	FORM ='UNFORMATTED', &
413	ACCESS ='DIRECT', &
414	RECL = 4096 )
415	IF( ios > 0 ) THEN
416	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obcsouth.dta file '
417	IF(lwp) WRITE(numout,*) '~~~~~~~'
418	nstop = nstop + 1
419	END IF
420	READ(inum,REC=1) clversion,clcom,irecl
421	CLOSE(inum)
422	IF(lwp) WRITE(numout,*)' '
423	IF(lwp) WRITE(numout,*)' opening obcsouth.dta with irecl=',irecl
424	OPEN(UNIT = inum, &
425	IOSTAT = ios, &
426	FILE ='obcsouth.dta', &
427	FORM ='UNFORMATTED', &
428	ACCESS ='DIRECT', &
429	RECL = irecl )
430	IF( ios > 0 ) THEN
431	IF(lwp) WRITE(numout,*) 'obc_dta: Pbm to OPEN the obcsouth.dta file '
432	IF(lwp) WRITE(numout,*) '~~~~~~~'
433	nstop = nstop + 1
434	END IF
435
436	! ... Read datafile and set temperature, salinity and normal velocity
437	! ... initialise the ssdta, tsdta arrays
438	! ... index 1 refer to before, 2 to after
439	DO jk = 1, jpkm1
440	irec = 2 + (jk -1)* jpf
441	READ(inum,REC=irec )((ssdta(jj,jk,1),ji=1,1),jj=jpisd, jpisf)
442	READ(inum,REC=irec+1)((tsdta(jj,jk,1),ji=1,1),jj=jpisd, jpisf)
443	READ(inum,REC=irec+2)((vsdta(jj,jk,1),ji=1,1),jj=jpisd, jpisf)
444	! ... do not interpolate: impose values of climatology
445	DO ji = 1, jpi
446	ii = ji -1 + nimpp
447	sfos(ji,jk) = ssdta(ii,jk,1)*tsmsk(ji,jk)
448	tfos(ji,jk) = tsdta(ii,jk,1)*tsmsk(ji,jk)
449	END DO
450	END DO
451	CLOSE(inum)
452
453	# if ! defined key_dynspg_fsc
454	! ... Rigid lid case: make sure vsdta is baroclinic velocity
455	! ... In rigid lid case vsdta needs to be the baroclinic component.
456	! ... substract the barotropic velocity component (zvsbtpe)
457	! ... from the total one (vsdta).
458
459	CALL obc_cli( vsdta, vclis, fs_njs0, fs_njs1, 1, jpi, nimpp )
460	# endif
461	! ... Set normal velocity (on njs0, njs1 <=> jpjsob)
462	DO jk = 1, jpkm1
463	DO ji = 1, jpi
464	ii = ji -1 + nimpp
465	vfos(ji,jk) = vsdta(ii,jk,1)*vsmsk(ji,jk)
466	END DO
467	END DO
468
469	END IF
470
471	ELSE
472
473	! 2. CALL not the first step ...
474	! do not read but interpolate between months.
475	! here no interpolation is done but the code is kept as a reminder
476	! ----------------------------------------------------------------------
477
478	IF( lpeastobc) THEN
479	DO jk = 1, jpkm1
480	DO jj = 1, jpj
481	ij = jj -1 + njmpp
482	sfoe(jj,jk) = sedta(ij,jk,1)*temsk(jj,jk)
483	tfoe(jj,jk) = tedta(ij,jk,1)*temsk(jj,jk)
484	ufoe(jj,jk) = uedta(ij,jk,1)*uemsk(jj,jk)
485	END DO
486	END DO
487	END IF
488
489	IF( lpwestobc) THEN
490	DO jk = 1, jpkm1
491	DO jj = 1, jpj
492	ij = jj -1 + njmpp
493	sfow(jj,jk) = swdta(ij,jk,1)*twmsk(jj,jk)
494	tfow(jj,jk) = twdta(ij,jk,1)*twmsk(jj,jk)
495	ufow(jj,jk) = uwdta(ij,jk,1)*uwmsk(jj,jk)
496	END DO
497	END DO
498	END IF
499
500	IF( lpnorthobc) THEN
501	DO jk = 1, jpkm1
502	DO ji = 1, jpi
503	ii = ji -1 + nimpp
504	sfon(ji,jk) = sndta(ii,jk,1)*tnmsk(ji,jk)
505	tfon(ji,jk) = tndta(ii,jk,1)*tnmsk(ji,jk)
506	vfon(ji,jk) = vndta(ii,jk,1)*vnmsk(ji,jk)
507	END DO
508	END DO
509	END IF
510
511	IF( lpsouthobc) THEN
512	DO jk = 1, jpkm1
513	DO ji = 1, jpi
514	ii = ji -1 + nimpp
515	sfos(ji,jk) = ssdta(ii,jk,1)*tsmsk(ji,jk)
516	tfos(ji,jk) = tsdta(ii,jk,1)*tsmsk(ji,jk)
517	vfos(ji,jk) = vsdta(ii,jk,1)*vsmsk(ji,jk)
518	END DO
519	END DO
520	END IF
521
522	END IF
523
524	END SUBROUTINE obc_dta
525
526	# if defined key_dynspg_fsc
527	!!-----------------------------------------------------------------------------
528	!! 'key_dynspg_fsc' free surface with constant volume
529	!!-----------------------------------------------------------------------------
530	SUBROUTINE obc_dta_psi ( kt ) ! Empty routine
531	!! * Arguments
532	INTEGER,INTENT(in) :: kt
533	WRITE(,) 'obc_dta_psi: You should not have seen this print! error?', kt
534	END SUBROUTINE obc_dta_psi
535	#else
536	!!-----------------------------------------------------------------------------
537	!! Default option Rigid-lid
538	!!-----------------------------------------------------------------------------
539
540	SUBROUTINE obc_dta_psi ( kt )
541	!!-----------------------------------------------------------------------------
542	!! * SUBROUTINE obc_dta_psi *
543	!!
544	!! ** Purpose :
545	!! Update the climatological streamfunction OBC at each time step.
546	!! Depends on the user's configuration. Here data are read only once
547	!! at the beginning of the run.
548	!!
549	!! ** Method :
550	!! 1. initialization
551	!! kbsfstart: number of time steps over which increase bsf
552	!! during initialization. This is provided for a smooth start
553	!! in cases where the transport is large (like on the Antarctic
554	!! continent). also note that when kbfstart=1, the transport
555	!! increases a lot in one time step and the precision usually
556	!! required for the solver may not be enough.
557	!! 2. set the time evolution of the climatological barotropic streamfunction
558	!! along the isolated coastlines ( gcbic(jnic) ).
559	!! 3. set the climatological barotropic streamfunction at the boundary.
560	!!
561	!! The last two steps are done only at first step (nit000) or if kt <= kbfstart
562	!!
563	!! History :
564	!! ! 97-08 (G. Madec, J.M. Molines)
565	!! 8.5 ! 02-10 (C. Talandier, A-M. Treguier) Free surface, F90
566	!!----------------------------------------------------------------------------
567	!! * Arguments
568	INTEGER, INTENT( in ) :: kt ! ocean time-step index
569
570	!! * Local declarations
571	INTEGER :: ji, jj, jnic, jip ! dummy loop indices
572	INTEGER :: inum = 11 ! temporary logical unit
573	INTEGER :: ip, ii, ij, iii, ijj
574	INTEGER :: kbsfstart
575	REAL(wp) :: zsver1, zsver2, zsver3, z2dtr, zcoef
576	!!----------------------------------------------------------------------------
577
578	! 1. initialisation
579	! -----------------
580
581	kbsfstart = 1
582	zsver1 = bsfic0(1)
583	zsver2 = zsver1
584	IF( kt <= kbsfstart ) THEN
585	zcoef = float(kt)/float(kbsfstart)
586	ELSE
587	zcoef = 1.
588	END IF
589	bsfic(1) = zsver1*zcoef
590	IF( lwp .AND. ( kt <= kbsfstart ) ) THEN
591	IF(lwp) WRITE(numout,*)' '
592	IF(lwp) WRITE(numout,*)'obcdta: spinup phase in obc_dta_psi routine'
593	IF(lwp) WRITE(numout,*)'~~~~~~ it=',kt,' OBC: spinup coef: ', &
594	zcoef, ' and transport: ',bsfic(1)
595	END IF
596
597	zsver2 = bsfic(1)-bsfic(2)
598	zsver3 = bsfic(2)
599
600	! 2. Right hand side of the barotropic elliptic equation (isolated coastlines)
601	! ----------------------------------------------------------------------------
602
603	IF( ( neuler == 0 ) .AND. ( kt == nit000 ) ) THEN
604	z2dtr = 1./rdt
605	ELSE
606	z2dtr = 1./2./rdt
607	END IF
608	! ... bsfb(ii,ij) should be constant but due to the Asselin filter it
609	! ... converges asymptotically towards bsfic(jnic)
610	! ... However, bsfb(ii,ij) is constant along the same coastlines
611	! ... ---> can be improved using an extra array for storing bsficb (before)
612	IF( nbobc > 1 ) THEN
613	DO jnic = 1,nbobc - 1
614	gcbic(jnic) = 0.
615	ip=mnic(0,jnic)
616	DO jip = 1,ip
617	ii = miic(jip,0,jnic)
618	ij = mjic(jip,0,jnic)
619	IF( ii >= nldi+ nimpp - 1 .AND. ii <= nlei+ nimpp - 1 .AND. &
620	ij >= nldj+ njmpp - 1 .AND. ij <= nlej+ njmpp - 1 ) THEN
621	iii=ii-nimpp+1
622	ijj=ij-njmpp+1
623	gcbic(jnic) = ( bsfic(jnic) - bsfb(iii,ijj) ) * z2dtr
624	END IF
625	END DO
626	END DO
627	END IF
628
629	IF( lk_mpp ) CALL mpp_isl( gcbic, 3 )
630
631	! 3. Update the climatological barotropic function at the boundary
632	! ----------------------------------------------------------------
633
634	IF( lpeastobc ) THEN
635
636	IF( kt == nit000 .OR. kt <= kbsfstart ) THEN
637	OPEN(inum,file='obceastbsf.dta')
638	READ(inum,*)
639	READ(inum,*)
640	READ(inum,*)
641	READ(inum,*)
642	READ(inum,*)
643	READ(inum,*) (bfoe(jj),jj=jpjed, jpjef)
644	CLOSE(inum)
645	END IF
646	DO jj=jpjed, jpjefm1
647	bfoe(jj)=bfoe(jj)*zcoef
648	END DO
649
650	END IF
651
652	IF( lpwestobc) THEN
653
654	IF( kt == nit000 .OR. kt <= kbsfstart ) then
655	OPEN(inum,file='obcwestbsf.dta')
656	READ(inum,*)
657	READ(inum,*)
658	READ(inum,*)
659	READ(inum,*)
660	READ(inum,*)
661	READ(inum,*) (bfow(jj),jj=jpjwd, jpjwf)
662	CLOSE(inum)
663	END IF
664	DO jj=jpjwd, jpjwfm1
665	bfow(jj)=bfow(jj)*zcoef
666	END DO
667
668	END IF
669
670	IF( lpsouthobc) THEN
671
672	IF( kt == nit000.OR.kt <= kbsfstart ) THEN
673	OPEN(inum,file='obcsouthbsf.dta')
674	READ(inum,*)
675	READ(inum,*)
676	READ(inum,*)
677	READ(inum,*)
678	READ(inum,*)
679	READ(inum,*) (bfos(jj),jj=jpisd, jpisf)
680	CLOSE(inum)
681	END IF
682	DO ji=jpisd, jpisfm1
683	bfos(ji)=bfos(ji)*zcoef
684	END DO
685
686	END IF
687
688	IF( lpnorthobc) THEN
689	IF( kt == nit000.OR.kt <= kbsfstart ) THEN
690	OPEN(inum,file='obcnorthbsf.dta')
691	READ(inum,*)
692	READ(inum,*)
693	READ(inum,*)
694	READ(inum,*)
695	READ(inum,*)
696	READ(inum,*) (bfon(jj),jj=jpind, jpinf)
697	CLOSE(inum)
698	END IF
699	DO ji=jpind, jpinfm1
700	bfon(ji)=bfon(ji)*zcoef
701	END DO
702
703	END IF
704
705	END SUBROUTINE obc_dta_psi
706
707	# endif
708
709	#else
710	!!------------------------------------------------------------------------------
711	!! default option: Dummy module NO Open Boundary Conditions
712	!!------------------------------------------------------------------------------
713	CONTAINS
714	SUBROUTINE obc_dta( kt ) ! Dummy routine
715	INTEGER, INTENT (in) :: kt
716	WRITE(,) 'obc_dta: You should not have seen this print! error?', kt
717	END SUBROUTINE obc_dta
718	#endif
719
720	!!==============================================================================
721	END MODULE obcdta

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

Download in other formats: