Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 2126)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 2128)
@@ -8,4 +8,6 @@
    !!             -   !  2007-07  (D. Storkey) add bdy_dta_bt
    !!            3.0  !  2008-04  (NEMO team)  add in the reference version
+   !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations 
+   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -25,4 +27,7 @@
    USE ioipsl
    USE in_out_manager  ! I/O logical units
+#if defined key_lim2
+   USE ice_2
+#endif
 
    IMPLICIT NONE
@@ -47,4 +52,10 @@
    REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta !: Arrays used for time interpolation of bdy data   
    REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            !: bdy data of ssh
+
+#if defined key_lim2
+   REAL(wp), DIMENSION(jpbdim,2)     ::   frld_bdydta          !: }
+   REAL(wp), DIMENSION(jpbdim,2)     ::   hicif_bdydta         !: } Arrays used for time interpolation of bdy data for ice variables
+   REAL(wp), DIMENSION(jpbdim,2)     ::   hsnif_bdydta         !: }
+#endif
 
    !!----------------------------------------------------------------------
@@ -196,5 +207,5 @@
                IF(lwp) WRITE(numout,*) 'offset: ',zoffset
                IF(lwp) WRITE(numout,*) 'totime: ',totime
-               IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr
+               IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr(1:ntimes_bdy)
 
                ! Check that there are not too many times in the file. 
@@ -227,17 +238,21 @@
                END IF
                !
-               IF    ( igrd == 1 ) THEN
-                 ntimes_bdyt = ntimes_bdy
-                 zoffsett = zoffset
-                 istept(:) = INT( zstepr(:) + zoffset )
-               ELSEIF(igrd == 2 ) THEN
-                 ntimes_bdyu = ntimes_bdy
-                 zoffsetu = zoffset
-                 istepu(:) = INT( zstepr(:) + zoffset )
-               ELSEIF(igrd == 3 ) THEN
-                 ntimes_bdyv = ntimes_bdy
-                 zoffsetv = zoffset
-                 istepv(:) = INT( zstepr(:) + zoffset )
-               ENDIF
+               SELECT CASE( igrd )
+                  CASE (1)
+                    ntimes_bdyt = ntimes_bdy
+                    zoffsett = zoffset
+                    istept(:) = INT( zstepr(:) + zoffset )
+                    numbdyt = inum
+                  CASE (2)
+                    ntimes_bdyu = ntimes_bdy
+                    zoffsetu = zoffset
+                    istepu(:) = INT( zstepr(:) + zoffset )
+                    numbdyu = inum
+                  CASE (3)
+                    ntimes_bdyv = ntimes_bdy
+                    zoffsetv = zoffset
+                    istepv(:) = INT( zstepr(:) + zoffset )
+                    numbdyv = inum
+               END SELECT
                !
             END DO                                         ! end loop over T, U & V grid 
@@ -291,6 +306,6 @@
             nbdy_b = it
             !
-            WRITE(numout,*) 'Time offset is ',zoffset
-            WRITE(numout,*) 'First record to read is ',nbdy_b
+            IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset
+            IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b
 
          ENDIF ! endif (nbdy_dta == 1)
@@ -334,4 +349,14 @@
             ENDIF
             !
+#if defined key_lim2
+            IF (ln_bdy_ice_frs) THEN
+               igrd = 1            ! T-points data
+               DO ib = 1, nblen(igrd)
+                  frld_bdy(ib)  =  frld(nbi(ib,igrd), nbj(ib,igrd))
+                  hicif_bdy(ib) = hicif(nbi(ib,igrd), nbj(ib,igrd))
+                  hsnif_bdy(ib) = hsnif(nbi(ib,igrd), nbj(ib,igrd))
+               END DO
+            ENDIF
+#endif
          ELSEIF( nbdy_dta == 1 ) THEN    ! Set first record in the climatological case:   
             !
@@ -357,5 +382,5 @@
                   nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
                ENDIF
-               WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd)
+               IF(lwp) WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd)
                ipi = nblendta(igrd)
                CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
@@ -372,5 +397,5 @@
                   nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
                ENDIF
-               WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd)
+               IF(lwp) WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd)
                ipi = nblendta(igrd)
                CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
@@ -390,5 +415,5 @@
                  nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar)
                ENDIF
-               WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd)
+               IF(lwp) WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd)
                ipi = nblendta(igrd)
                CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
@@ -404,5 +429,5 @@
                  nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar)
                ENDIF
-               WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd)
+               IF(lwp) WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd)
                ipi = nblendta(igrd)
                CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
@@ -414,4 +439,32 @@
             ENDIF ! ln_bdy_dyn_frs
 
+#if defined key_lim2
+            IF(ln_bdy_ice_frs) THEN
+              !
+              igrd=1                                              ! leads fraction
+              IF(lwp) WRITE(numout,*) 'Dim size for ildsconc is ',nblendta(igrd)
+              ipi=nblendta(igrd)
+              CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a )
+              DO ib=1, nblen(igrd)
+                frld_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
+              END DO
+              !
+              igrd=1                                              ! ice thickness
+              IF(lwp) WRITE(numout,*) 'Dim size for iicethic is ',nblendta(igrd)
+              ipi=nblendta(igrd)
+              CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a )
+              DO ib=1, nblen(igrd)
+                hicif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
+              END DO
+              !
+              igrd=1                                              ! snow thickness
+              IF(lwp) WRITE(numout,*) 'Dim size for isnowthi is ',nblendta(igrd)
+              ipi=nblendta(igrd)
+              CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a )
+              DO ib=1, nblen(igrd)
+                hsnif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
+              END DO
+            ENDIF ! just if ln_bdy_ice_frs is set
+#endif
 
             IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN
@@ -427,4 +480,11 @@
                  vbdydta(:,:,1) = vbdydta(:,:,2)
                ENDIF
+#if defined key_lim2
+             IF( ln_bdy_ice_frs ) THEN
+               frld_bdydta (:,1) =  frld_bdydta(:,2)
+               hicif_bdydta(:,1) = hicif_bdydta(:,2)
+               hsnif_bdydta(:,1) = hsnif_bdydta(:,2)
+             ENDIF
+#endif
             END IF
 
@@ -441,4 +501,11 @@
               vbdy  (:,:) = vbdydta  (:,:,2)
             ENDIF
+#if defined key_lim2
+          IF(ln_bdy_ice_frs) THEN
+            frld_bdy (:) = frld_bdydta (:,2)
+            hicif_bdy(:) = hicif_bdydta(:,2)
+            hsnif_bdy(:) = hsnif_bdydta(:,2)
+          ENDIF
+#endif
 
             IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )
@@ -483,4 +550,11 @@
                 vbdydta(:,:,1) =  vbdydta(:,:,2)
               ENDIF
+#if defined key_lim2
+              IF(ln_bdy_ice_frs) THEN
+                frld_bdydta (:,1) =  frld_bdydta (:,2)
+                hicif_bdydta(:,1) =  hicif_bdydta(:,2)
+                hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
+              ENDIF
+#endif
             END IF ! nbdy_a < ntimes_bdy
 
@@ -497,5 +571,11 @@
              vbdydta(:,:,1) =  vbdydta(:,:,2)
            ENDIF
- 
+#if defined key_lim2
+           IF(ln_bdy_ice_frs) THEN
+             frld_bdydta (:,1) =  frld_bdydta (:,2)
+             hicif_bdydta(:,1) =  hicif_bdydta(:,2)
+             hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
+           ENDIF
+#endif 
            ! read another set
            ipj  = 1
@@ -543,5 +623,30 @@
               END DO
            ENDIF ! ln_bdy_dyn_frs
-
+           !
+#if defined key_lim2
+           IF(ln_bdy_ice_frs) THEN
+             !
+             igrd = 1                                    ! ice concentration
+             ipi=nblendta(igrd)
+             CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a )
+             DO ib=1, nblen(igrd)
+               frld_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
+             END DO
+             !
+             igrd=1                                      ! ice thickness
+             ipi=nblendta(igrd)
+             CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a )
+             DO ib=1, nblen(igrd)
+               hicif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
+             END DO
+             !
+             igrd=1                                      ! snow thickness
+             ipi=nblendta(igrd)
+             CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a )
+             DO ib=1, nblen(igrd)
+               hsnif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
+             END DO
+           ENDIF ! ln_bdy_ice_frs
+#endif
            !
            IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b
@@ -559,6 +664,8 @@
        ! ********************
        ! 
-       IF( ln_bdy_clim ) THEN   ;   zxy = REAL( nday                  , wp ) / REAL( nmonth_len(nbdy_b), wp ) + 0.5 - i15
-       ELSE                     ;   zxy = REAL( istep(nbdy_b) - itimer, wp ) / REAL( istep(nbdy_b) - istep(nbdy_a), wp )
+       IF( ln_bdy_clim ) THEN   ;   zxy = REAL( nday                   ) / REAL( nmonth_len(nbdy_b) ) + 0.5 - i15
+       ELSEIF( istep(nbdy_b) == istep(nbdy_a) ) THEN 
+                                    zxy = 0.0_wp
+       ELSE                     ;   zxy = REAL( istep(nbdy_b) - itimer ) / REAL( istep(nbdy_b) - istep(nbdy_a) )
        END IF
 
@@ -589,4 +696,15 @@
           ENDIF
 
+#if defined key_lim2
+          IF(ln_bdy_ice_frs) THEN
+            igrd=1
+            DO ib=1, nblen(igrd)
+               frld_bdy(ib) = zxy *  frld_bdydta(ib,2) + (1.-zxy) *  frld_bdydta(ib,1)
+              hicif_bdy(ib) = zxy * hicif_bdydta(ib,2) + (1.-zxy) * hicif_bdydta(ib,1)
+              hsnif_bdy(ib) = zxy * hsnif_bdydta(ib,2) + (1.-zxy) * hsnif_bdydta(ib,1)
+            END DO
+          ENDIF ! just if ln_bdy_ice_frs is true
+#endif
+
       END IF                       !end if ((nbdy_dta==1).AND.(ntimes_bdy>1))
     
@@ -606,5 +724,5 @@
 
 
-   SUBROUTINE bdy_dta_bt( kt, jit )
+   SUBROUTINE bdy_dta_bt( kt, jit, icycl )
       !!---------------------------------------------------------------------------
       !!                      ***  SUBROUTINE bdy_dta_bt  ***
@@ -620,4 +738,5 @@
       INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
       INTEGER, INTENT( in ) ::   jit         ! barotropic time step index
+      INTEGER, INTENT( in ) ::   icycl       ! number of cycles need for final file close
       !                                      ! (for timesplitting option, otherwise zero)
       !!
@@ -639,5 +758,5 @@
       REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files
       REAL(wp), DIMENSION(jpbdta,1)     ::   zdta               ! temporary array for data fields
-      CHARACTER(LEN=80), DIMENSION(3)   ::   clfile
+      CHARACTER(LEN=80), DIMENSION(6)   ::   clfile
       CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate
       !!---------------------------------------------------------------------------
@@ -688,5 +807,5 @@
 
       !                                                !-------------------!
-      IF( kt == nit000 ) THEN                          !  First call only  !
+      IF( kt == nit000 .and. jit ==2 ) THEN            !  First call only  !
          !                                             !-------------------!
          istep_bt(:) = 0
@@ -712,20 +831,20 @@
                                                      ! necessary time dumps in file are included
 
-          clfile(1) = filbdy_data_bt_T
-          clfile(2) = filbdy_data_bt_U
-          clfile(3) = filbdy_data_bt_V
-
-          DO igrd = 1,3
+          clfile(4) = filbdy_data_bt_T
+          clfile(5) = filbdy_data_bt_U
+          clfile(6) = filbdy_data_bt_V
+
+          DO igrd = 4,6
 
             CALL iom_open( clfile(igrd), inum )
-            CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 
+            CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy_bt, cdunits=clunits ) 
 
             SELECT CASE( igrd )
-               CASE (1) 
-                  numbdyt = inum
-               CASE (2) 
-                  numbdyu = inum
-               CASE (3) 
-                  numbdyv = inum
+               CASE (4) 
+                  numbdyt_bt = inum
+               CASE (5) 
+                  numbdyu_bt = inum
+               CASE (6) 
+                  numbdyv_bt = inum
             END SELECT
 
@@ -757,5 +876,5 @@
 
             ! Check that time array increases (or interp will fail):
-            DO it = 2, ntimes_bdy
+            DO it = 2, ntimes_bdy_bt
                IF ( zstepr(it-1) >= zstepr(it) ) THEN
                   CALL ctl_stop('Time array in unstructured boundary data file', &
@@ -778,38 +897,30 @@
                ! The same applies to the last time level: see setting of lect below.
 
-               IF ( ntimes_bdy == 1 ) CALL ctl_stop( &
+               IF ( ntimes_bdy_bt == 1 ) CALL ctl_stop( &
                     'There is only one time dump in data files', &
                     'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
 
                zinterval_s = zstepr(2) - zstepr(1)
-               zinterval_e = zstepr(ntimes_bdy) - zstepr(ntimes_bdy-1)
-
-               IF ( zstepr(1) - zinterval_s / 2.0 > 0 ) THEN	          
-                  IF(lwp) WRITE(numout,*) 'First bdy time relative to nit000:', zstepr(1)
-                  IF(lwp) WRITE(numout,*) 'Interval between first two times: ', zinterval_s
-                  CALL ctl_stop( 'First data time is after start of run', & 
-                       'by more than half a meaning period', &
-                       'Check file: ' // TRIM(clfile(igrd)) )
+               zinterval_e = zstepr(ntimes_bdy_bt) - zstepr(ntimes_bdy_bt-1)
+
+               IF( zstepr(1) + zoffset > 0 ) THEN
+                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
+                     CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 )
                END IF
-
-               IF ( zstepr(ntimes_bdy) + zinterval_e / 2.0 < totime ) THEN
-                  IF(lwp) WRITE(numout,*) 'Last bdy time relative to nit000:', zstepr(ntimes_bdy)
-                  IF(lwp) WRITE(numout,*) 'Interval between last two times: ', zinterval_e
-                  CALL ctl_stop( 'Last data time is before end of run', & 
-                       'by more than half a meaning period', &
-                       'Check file: ' // TRIM(clfile(igrd))  )
+               IF( zstepr(ntimes_bdy_bt) + zoffset < totime ) THEN
+                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
+                     CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 )
                END IF
-
             END IF ! .NOT. ln_bdy_clim
 
-            IF ( igrd .EQ. 1) THEN
+            IF ( igrd .EQ. 4) THEN
               ntimes_bdyt = ntimes_bdy_bt
               zoffsett = zoffset
               istept(:) = INT( zstepr(:) + zoffset )
-            ELSE IF (igrd .EQ. 2) THEN
+            ELSE IF (igrd .EQ. 5) THEN
               ntimes_bdyu = ntimes_bdy_bt
               zoffsetu = zoffset
               istepu(:) = INT( zstepr(:) + zoffset )
-            ELSE IF (igrd .EQ. 3) THEN
+            ELSE IF (igrd .EQ. 6) THEN
               ntimes_bdyv = ntimes_bdy_bt
               zoffsetv = zoffset
@@ -865,6 +976,6 @@
           nbdy_b_bt = it
 
-          WRITE(numout,*) 'Time offset is ',zoffset
-          WRITE(numout,*) 'First record to read is ',nbdy_b_bt
+          IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset
+          IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b_bt
 
         ENDIF ! endif (nbdy_dta == 1)
@@ -875,15 +986,15 @@
         IF ( nbdy_dta == 0) THEN
           ! boundary data arrays are filled with initial conditions
-          igrd = 2            ! U-points data 
+          igrd = 5            ! U-points data 
           DO ib = 1, nblen(igrd)              
             ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1)
           END DO
 
-          igrd = 3            ! V-points data 
+          igrd = 6            ! V-points data 
           DO ib = 1, nblen(igrd)              
             vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1)
           END DO
 
-          igrd = 1            ! T-points data 
+          igrd = 4            ! T-points data 
           DO ib = 1, nblen(igrd)              
             sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd))
@@ -910,9 +1021,9 @@
          ! Read first record:
           ipj=1
-          igrd=1
+          igrd=4
           ipi=nblendta(igrd)
 
           ! ssh
-          igrd=1
+          igrd=4
           IF ( nblendta(igrd) .le. 0 ) THEN 
             idvar = iom_varid( numbdyt_bt,'sossheig' )
@@ -929,5 +1040,5 @@
  
           ! u-velocity
-          igrd=2
+          igrd=5
           IF ( nblendta(igrd) .le. 0 ) THEN 
             idvar = iom_varid( numbdyu_bt,'vobtcrtx' )
@@ -944,5 +1055,5 @@
 
           ! v-velocity
-          igrd=3
+          igrd=6
           IF ( nblendta(igrd) .le. 0 ) THEN 
             idvar = iom_varid( numbdyv_bt,'vobtcrty' )
@@ -1010,10 +1121,10 @@
           ipj=1
           ipk=jpk
-          igrd=1
+          igrd=4
           ipi=nblendta(igrd)
 
           
           ! ssh
-          igrd=1
+          igrd=4
           ipi=nblendta(igrd)
 
@@ -1025,5 +1136,5 @@
 
           ! u-velocity
-          igrd=2
+          igrd=5
           ipi=nblendta(igrd)
 
@@ -1035,5 +1146,5 @@
 
           ! v-velocity
-          igrd=3
+          igrd=6
           ipi=nblendta(igrd)
 
@@ -1064,5 +1175,5 @@
         END IF
 
-          igrd=1
+          igrd=4
           DO ib=1, nblen(igrd)
             sshbdy(ib) = zxy      * sshbdydta(ib,2) + &
@@ -1070,5 +1181,5 @@
           END DO
 
-          igrd=2
+          igrd=5
           DO ib=1, nblen(igrd)
             ubtbdy(ib) = zxy      * ubtbdydta(ib,2) + &
@@ -1076,5 +1187,5 @@
           END DO
 
-          igrd=3
+          igrd=6
           DO ib=1, nblen(igrd)
             vbtbdy(ib) = zxy      * vbtbdydta(ib,2) + &
@@ -1090,5 +1201,5 @@
 
       ! Closing of the 3 files
-      IF( kt == nitend ) THEN
+      IF( kt == nitend   .and. jit == icycl ) THEN
           CALL iom_close( numbdyt_bt )
           CALL iom_close( numbdyu_bt )
@@ -1109,5 +1220,5 @@
       WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dta
-   SUBROUTINE bdy_dta_bt( kt, kit )      ! Empty routine
+   SUBROUTINE bdy_dta_bt( kt, kit, icycle )      ! Empty routine
       WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, kit
    END SUBROUTINE bdy_dta_bt
