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

Last change on this file since 78 was 78, checked in by opalod, 20 years ago

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