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

source: trunk/NEMO/OPA_SRC/OBC/obcdta.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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