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

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

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