Changeset 8813 for branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC

Timestamp:

2017-11-24T17:56:51+01:00 (6 years ago)

Author:

gm

Message:

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8787

Location:

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC

Files:

• BDY/bdydta.F90 (modified) (45 diffs)
• SBC/sbc_ice.F90 (modified) (6 diffs)
• SBC/sbcblk.F90 (modified) (18 diffs)
• SBC/sbcblk_algo_coare.F90 (modified) (1 diff)
• SBC/sbccpl.F90 (modified) (17 diffs)
• SBC/sbcice_cice.F90 (modified) (1 diff)
• SBC/sbcisf.F90 (modified) (4 diffs)

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

@@ -51 +51 @@
 
 #if defined key_lim3
-   LOGICAL :: ll_bdylim3                  ! determine whether ice input is 1cat (F) or Xcat (T) type
-   INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure
+   INTEGER ::   nice_cat                      ! number of categories in the input file
+   INTEGER ::   jfld_hti, jfld_hts, jfld_ai   ! indices of ice thickness, snow thickness and concentration in bf structure
 #endif
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -80 +80 @@
       !                                               ! etc.
       !
-      INTEGER ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices
+      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl  ! dummy loop indices
+      INTEGER ::  ii, ij, ik, igrd                  ! local integers
       INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
       INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
@@ -95 +96 @@
          !-----------------------------
          
-         DO ib_bdy = 1, nb_bdy
+         DO jbdy = 1, nb_bdy
             !
-            nblen => idx_bdy(ib_bdy)%nblen
-            nblenrim => idx_bdy(ib_bdy)%nblenrim
-            dta => dta_bdy(ib_bdy)
-
-            IF( nn_dyn2d_dta(ib_bdy) == 0 ) THEN 
+            nblen    => idx_bdy(jbdy)%nblen
+            nblenrim => idx_bdy(jbdy)%nblenrim
+            dta      => dta_bdy(jbdy)
+            !
+            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_ssh ) THEN 
                   igrd = 1
                   DO ib = 1, ilen1(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_u2d ) THEN 
                   igrd = 2
                   DO ib = 1, ilen1(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
+                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_v2d ) THEN 
                   igrd = 3
                   DO ib = 1, ilen1(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
+                     ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
                   END DO 
-               END IF
-            ENDIF
-
-            IF( nn_dyn3d_dta(ib_bdy) == 0 ) THEN 
+               ENDIF
+            ENDIF
+            !
+            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_u3d ) THEN 
@@ -135 +136 @@
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
                      END DO
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_v3d ) THEN 
                   igrd = 3 
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
                         END DO
                   END DO 
-               END IF
-            ENDIF
-
-            IF( nn_tra_dta(ib_bdy) == 0 ) THEN 
+               ENDIF
+            ENDIF
+
+            IF( nn_tra_dta(jbdy) == 0 ) THEN 
                ilen1(:) = nblen(:)
                IF( dta%ll_tem ) THEN
@@ -159 +160 @@
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
                      END DO
                   END DO 
-               END IF
+               ENDIF
                IF( dta%ll_sal ) THEN
                   igrd = 1 
                   DO ib = 1, ilen1(igrd)
                      DO ik = 1, jpkm1
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
                      END DO
                   END DO 
-               END IF
-            ENDIF
-
-#if defined key_lim3
-            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN 
+               ENDIF
+            ENDIF
+
+#if defined key_lim3
+            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN    ! set ice to initial values
                ilen1(:) = nblen(:)
                IF( dta%ll_a_i ) THEN
@@ -184 +185 @@
                   DO jl = 1, jpl
                      DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -194 +195 @@
                   DO jl = 1, jpl
                      DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -204 +205 @@
                   DO jl = 1, jpl
                      DO ib = 1, ilen1(igrd)
-                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
+                        ii = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij = idx_bdy(jbdy)%nbj(ib,igrd)
+                        dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1) 
                      END DO
                   END DO
@@ -212 +213 @@
             ENDIF
 #endif
-         END DO ! ib_bdy
+         END DO ! jbdy
          !
       ENDIF ! kt == nit000
@@ -220 +221 @@
      
       jstart = 1
-      DO ib_bdy = 1, nb_bdy   
-         dta => dta_bdy(ib_bdy)
-         IF( nn_dta(ib_bdy) == 1 ) THEN ! skip this bit if no external data required
+      DO jbdy = 1, nb_bdy   
+         dta => dta_bdy(jbdy)
+         IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required
       
             IF( PRESENT(jit) ) THEN
                ! Update barotropic boundary conditions only
                ! jit is optional argument for fld_read and bdytide_update
-               IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN
-                  IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+               IF( cn_dyn2d(jbdy) /= 'none' ) THEN
+                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
                      IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
                      IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
                      IF( dta%ll_u3d ) dta%v2d(:) = 0._wp
                   ENDIF
-                  IF (cn_tra(ib_bdy) /= 'runoff') THEN
-                     IF( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 ) THEN
+                  IF (cn_tra(jbdy) /= 'runoff') THEN
+                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN
 
                         jend = jstart + dta%nread(2) - 1
-                        IF( ln_full_vel_array(ib_bdy) ) THEN
+                        IF( ln_full_vel_array(jbdy) ) THEN
                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy)  )
+                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(jbdy)  )
                         ELSE
                            CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
@@ -246 +247 @@
 
                         ! If full velocities in boundary data then extract barotropic velocities from 3D fields
-                        IF( ln_full_vel_array(ib_bdy) .AND.                                             &
-                          &    ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR.  &
-                          &      nn_dyn3d_dta(ib_bdy) == 1 ) )THEN
+                        IF( ln_full_vel_array(jbdy) .AND.                                             &
+                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  &
+                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN
 
                            igrd = 2                      ! zonal velocity
                            dta%u2d(:) = 0._wp
-                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                               DO ik = 1, jpkm1
                                  dta%u2d(ib) = dta%u2d(ib) &
@@ -263 +264 @@
                            igrd = 3                      ! meridional velocity
                            dta%v2d(:) = 0._wp
-                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                               DO ik = 1, jpkm1
                                  dta%v2d(ib) = dta%v2d(ib) &
@@ -274 +275 @@
                         ENDIF                    
                      ENDIF
-                     IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta, td=tides(ib_bdy),   & 
+                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   & 
                           &                 jit=jit, time_offset=time_offset )
                      ENDIF
@@ -281 +282 @@
                ENDIF
             ELSE
-               IF (cn_tra(ib_bdy) == 'runoff') then      ! runoff condition
-                  jend = nb_bdy_fld(ib_bdy)
+               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition
+                  jend = nb_bdy_fld(jbdy)
                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
                                & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
                   !
                   igrd = 2                      ! zonal velocity
-                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
                   END DO
                   !
                   igrd = 3                      ! meridional velocity
-                  DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                     ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                     ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
                   END DO
                ELSE
-                  IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
                      IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
                      IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
@@ -308 +309 @@
                      jend = jstart + dta%nread(1) - 1
                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
-                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) )
+                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(jbdy) )
                   ENDIF
                   ! If full velocities in boundary data then split into barotropic and baroclinic data
-                  IF( ln_full_vel_array(ib_bdy) .and.                                             &
-                    & ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR. &
-                    &   nn_dyn3d_dta(ib_bdy) == 1 ) ) THEN
+                  IF( ln_full_vel_array(jbdy) .and.                                             &
+                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. &
+                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN
                      igrd = 2                      ! zonal velocity
                      dta%u2d(:) = 0._wp
-                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                         DO ik = 1, jpkm1
                            dta%u2d(ib) = dta%u2d(ib) &
@@ -330 +331 @@
                      igrd = 3                      ! meridional velocity
                      dta%v2d(:) = 0._wp
-                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                         DO ik = 1, jpkm1
                            dta%v2d(ib) = dta%v2d(ib) &
@@ -346 +347 @@
                ENDIF
 #if defined key_lim3
-               IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is 1cat)
-                CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
-                                  & dta_bdy(ib_bdy)%h_i,     dta_bdy(ib_bdy)%h_s,     dta_bdy(ib_bdy)%a_i     )
+               IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1 ) THEN
+                  IF( nice_cat == 1 ) THEN ! case input cat = 1
+                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
+                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
+                  ELSEIF( nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl
+                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), &
+                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
+                  ENDIF
                ENDIF
 #endif
             ENDIF
             jstart = jstart + dta%nread(1)
-         END IF ! nn_dta(ib_bdy) = 1
-      END DO  ! ib_bdy
+         ENDIF    ! nn_dta(jbdy) = 1
+      END DO  ! jbdy
 
       IF ( ln_tide ) THEN
          IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                           
-            DO ib_bdy = 1, nb_bdy    ! Tidal component added in ts loop
-               IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
-                  nblen => idx_bdy(ib_bdy)%nblen
-                  nblenrim => idx_bdy(ib_bdy)%nblenrim
-                  IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
-                  IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen1(1)) = dta_bdy(ib_bdy)%ssh(1:ilen1(1))
-                  IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen1(2)) = dta_bdy(ib_bdy)%u2d(1:ilen1(2))
-                  IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen1(3)) = dta_bdy(ib_bdy)%v2d(1:ilen1(3))
+            DO jbdy = 1, nb_bdy    ! Tidal component added in ts loop
+               IF ( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN
+                  nblen => idx_bdy(jbdy)%nblen
+                  nblenrim => idx_bdy(jbdy)%nblenrim
+                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
+                  IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
+                  IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
+                  IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
                ENDIF
             END DO
@@ -375 +381 @@
 
       IF ( ln_apr_obc ) THEN
-         DO ib_bdy = 1, nb_bdy
-            IF (cn_tra(ib_bdy) /= 'runoff')THEN
+         DO jbdy = 1, nb_bdy
+            IF (cn_tra(jbdy) /= 'runoff')THEN
                igrd = 1                      ! meridional velocity
-               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij)
+               DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd)
+                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                  dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij)
                END DO
             ENDIF
@@ -402 +408 @@
       !!                
       !!----------------------------------------------------------------------
-      INTEGER ::   ib_bdy, jfld, jstart, jend, ierror, ios     ! Local integers
+      INTEGER ::   jbdy, jfld, jstart, jend, ierror, ios     ! Local integers
       !
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
@@ -408 +414 @@
       CHARACTER(len = 256)::   clname                           ! temporary file name
       LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
-                                                                ! =F => baroclinic velocities in 3D boundary data
+      !                                                         ! =F => baroclinic velocities in 3D boundary data
       INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays
@@ -416 +422 @@
       TYPE(OBC_DATA), POINTER                ::   dta           ! short cut
 #if defined key_lim3
-      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
+      INTEGER               ::   kndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat)
+      INTEGER, DIMENSION(4) ::   kdimsz   ! size   of dimensions
       INTEGER               ::   inum,id1 ! local integer
 #endif
@@ -440 +447 @@
 
       ! Set nn_dta
-      DO ib_bdy = 1, nb_bdy
-         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       &
-                               ,nn_dyn3d_dta(ib_bdy)       &
-                               ,nn_tra_dta(ib_bdy)         &
-#if defined key_lim3
-                              ,nn_ice_lim_dta(ib_bdy)    &
+      DO jbdy = 1, nb_bdy
+         nn_dta(jbdy) = MAX(   nn_dyn2d_dta  (jbdy)    &
+            &                , nn_dyn3d_dta  (jbdy)    &
+            &                , nn_tra_dta    (jbdy)    &
+#if defined key_lim3
+            &                , nn_ice_lim_dta(jbdy)    &
 #endif
                               )
-         IF(nn_dta(ib_bdy) > 1) nn_dta(ib_bdy) = 1
+         IF(nn_dta(jbdy) > 1)   nn_dta(jbdy) = 1
       END DO
 
@@ -455 +462 @@
       ALLOCATE( nb_bdy_fld(nb_bdy) )
       nb_bdy_fld(:) = 0
-      DO ib_bdy = 1, nb_bdy         
-         IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
-         ENDIF
-         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) == 1 ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
-         ENDIF
-         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) == 1  ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
-         ENDIF
-#if defined key_lim3
-         IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) == 1  ) THEN
-            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
+      DO jbdy = 1, nb_bdy         
+         IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
+         ENDIF
+         IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
+         ENDIF
+         IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1  ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2
+         ENDIF
+#if defined key_lim3
+         IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1  ) THEN
+            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3
          ENDIF
 #endif               
-         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy)
+         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy)
       END DO            
 
@@ -496 +503 @@
       REWIND(numnam_cfg)
       jfld = 0 
-      DO ib_bdy = 1, nb_bdy         
-         IF( nn_dta(ib_bdy) == 1 ) THEN
+      DO jbdy = 1, nb_bdy         
+         IF( nn_dta(jbdy) == 1 ) THEN
             READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp )
@@ -505 +512 @@
             IF(lwm) WRITE( numond, nambdy_dta )
 
-            cn_dir_array(ib_bdy) = cn_dir
-            ln_full_vel_array(ib_bdy) = ln_full_vel
-
-            nblen => idx_bdy(ib_bdy)%nblen
-            nblenrim => idx_bdy(ib_bdy)%nblenrim
-            dta => dta_bdy(ib_bdy)
+            cn_dir_array(jbdy) = cn_dir
+            ln_full_vel_array(jbdy) = ln_full_vel
+
+            nblen => idx_bdy(jbdy)%nblen
+            nblenrim => idx_bdy(jbdy)%nblenrim
+            dta => dta_bdy(jbdy)
             dta%nread(2) = 0
 
             ! Only read in necessary fields for this set.
             ! Important that barotropic variables come first.
-            IF( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN 
+            IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 
 
                IF( dta%ll_ssh ) THEN 
@@ -521 +528 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -528 +535 @@
                ENDIF
 
-               IF( dta%ll_u2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN
+               IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in u2d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u2d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 2
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -539 +546 @@
                ENDIF
 
-               IF( dta%ll_v2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN
+               IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in v2d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_v2d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 3
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -554 +561 @@
             ! read 3D velocities if baroclinic velocities require OR if
             ! barotropic velocities required and ln_full_vel set to .true.
-            IF( nn_dyn3d_dta(ib_bdy) == 1 .OR. &
-           &  ( ln_full_vel_array(ib_bdy) .AND. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN
-
-               IF( dta%ll_u3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
+            IF( nn_dyn3d_dta(jbdy) == 1 .OR. &
+           &  ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
+
+               IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in u3d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u3d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 2
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = jpk
-                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
-               ENDIF
-
-               IF( dta%ll_v3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
+                  IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1
+               ENDIF
+
+               IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
                   if(lwp) write(numout,*) '++++++ reading in v3d field'
                   jfld = jfld + 1
                   blf_i(jfld) = bn_v3d
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 3
                   ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = jpk
-                  IF( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
+                  IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1
                ENDIF
 
@@ -582 +589 @@
 
             ! temperature and salinity
-            IF( nn_tra_dta(ib_bdy) == 1 ) THEN
+            IF( nn_tra_dta(jbdy) == 1 ) THEN
 
                IF( dta%ll_tem ) THEN
@@ -588 +595 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_tem
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -598 +605 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_sal
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld) = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
@@ -608 +615 @@
 #if defined key_lim3
             ! sea ice
-            IF( nn_ice_lim_dta(ib_bdy) == 1 ) THEN
+            IF( nn_ice_lim_dta(jbdy) == 1 ) THEN
                ! Test for types of ice input (1cat or Xcat) 
                ! Build file name to find dimensions 
@@ -622 +629 @@
                !
                CALL iom_open  ( clname, inum )
-               id1 = iom_varid( inum, bn_a_i%clvar, kndims=zndims, ldstop = .FALSE. )
+               id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. )
                CALL iom_close ( inum )
 
-                IF ( zndims == 4 ) THEN
-                 ll_bdylim3 = .TRUE.   ! Xcat input
+                IF ( kndims == 4 ) THEN
+                 nice_cat = kdimsz(4)   ! Xcat input
                ELSE
-                 ll_bdylim3 = .FALSE.  ! 1cat input      
+                 nice_cat = 1           ! 1cat input      
                ENDIF
                ! End test
@@ -635 +642 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_a_i
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld)  = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
-                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+                  ilen3(jfld) = nice_cat
                ENDIF
 
@@ -644 +651 @@
                   jfld = jfld + 1
                   blf_i(jfld) = bn_h_i
-                  ibdy(jfld) = ib_bdy
+                  ibdy(jfld)  = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
-                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+                  ilen3(jfld) = nice_cat
                ENDIF
 
                IF( dta%ll_h_s ) THEN
                   jfld = jfld + 1
-                   blf_i(jfld) = bn_h_s
-                  ibdy(jfld) = ib_bdy
+                  blf_i(jfld) = bn_h_s
+                  ibdy(jfld)  = jbdy
                   igrid(jfld) = 1
                   ilen1(jfld) = nblen(igrid(jfld))
-                  IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF
+                  ilen3(jfld) = nice_cat
                ENDIF
 
@@ -663 +670 @@
             ! Recalculate field counts
             !-------------------------
-            IF( ib_bdy == 1 ) THEN 
+            IF( jbdy == 1 ) THEN 
                nb_bdy_fld_sum = 0
-               nb_bdy_fld(ib_bdy) = jfld
+               nb_bdy_fld(jbdy) = jfld
                nb_bdy_fld_sum     = jfld              
             ELSE
-               nb_bdy_fld(ib_bdy) = jfld - nb_bdy_fld_sum
-               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(ib_bdy)
-            ENDIF
-
-            dta%nread(1) = nb_bdy_fld(ib_bdy)
+               nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum
+               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy)
+            ENDIF
+
+            dta%nread(1) = nb_bdy_fld(jbdy)
 
          ENDIF ! nn_dta == 1
-      ENDDO ! ib_bdy
+      ENDDO ! jbdy
 
       DO jfld = 1, nb_bdy_fld_sum
@@ -687 +694 @@
       !-------------------------------------
       jstart = 1
-      DO ib_bdy = 1, nb_bdy
-         jend = jstart - 1 + nb_bdy_fld(ib_bdy) 
-         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   &
+      DO jbdy = 1, nb_bdy
+         jend = jstart - 1 + nb_bdy_fld(jbdy) 
+         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta',   &
          &              'open boundary conditions', 'nambdy_dta' )
          jstart = jend + 1
@@ -700 +707 @@
 
       jfld = 0
-      DO ib_bdy=1, nb_bdy
-
-         nblen => idx_bdy(ib_bdy)%nblen
-         dta => dta_bdy(ib_bdy)
+      DO jbdy=1, nb_bdy
+
+         nblen => idx_bdy(jbdy)%nblen
+         dta => dta_bdy(jbdy)
 
          if(lwp) then
@@ -715 +722 @@
          endif
 
-         IF ( nn_dyn2d_dta(ib_bdy) == 0 .or. nn_dyn2d_dta(ib_bdy) == 2 ) THEN
+         IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN
             if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space'
             IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) )
@@ -721 +728 @@
             IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) )
          ENDIF
-         IF ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN
+         IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN
             IF( dta%ll_ssh ) THEN
                if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow'
@@ -728 +735 @@
             ENDIF
             IF ( dta%ll_u2d ) THEN
-               IF ( ln_full_vel_array(ib_bdy) ) THEN
+               IF ( ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ dta%u2d allocated space'
                   ALLOCATE( dta%u2d(nblen(2)) )
@@ -738 +745 @@
             ENDIF
             IF ( dta%ll_v2d ) THEN
-               IF ( ln_full_vel_array(ib_bdy) ) THEN
+               IF ( ln_full_vel_array(jbdy) ) THEN
                   if(lwp) write(numout,*) '++++++ dta%v2d allocated space'
                   ALLOCATE( dta%v2d(nblen(3)) )
@@ -749 +756 @@
          ENDIF
 
-         IF ( nn_dyn3d_dta(ib_bdy) == 0 ) THEN
+         IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN
             if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space'
-            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) )
-            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) )
-         ENDIF
-         IF ( nn_dyn3d_dta(ib_bdy) == 1 .or. &
-           &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN
-            IF ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN
+            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) )
+            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) )
+         ENDIF
+         IF ( nn_dyn3d_dta(jbdy) == 1 .or. &
+           &  ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN
+            IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN
                if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:)
-            ENDIF
-            IF ( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN
+               dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:)
+            ENDIF
+            IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN
                if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:)
-            ENDIF
-         ENDIF
-
-         IF( nn_tra_dta(ib_bdy) == 0 ) THEN
+               dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:)
+            ENDIF
+         ENDIF
+
+         IF( nn_tra_dta(jbdy) == 0 ) THEN
             if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space'
-            IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) )
-            IF( dta%ll_sal ) ALLOCATE( dta_bdy(ib_bdy)%sal(nblen(1),jpk) )
+            IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) )
+            IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) )
          ELSE
             IF( dta%ll_tem ) THEN
                if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:)
+               dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:)
             ENDIF
             IF( dta%ll_sal ) THEN 
                if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow'
                jfld = jfld + 1
-               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:)
-            ENDIF
-         ENDIF
-
-#if defined key_lim3
-         IF (cn_ice_lim(ib_bdy) /= 'none') THEN
-            IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN
-               ALLOCATE( dta_bdy(ib_bdy)%a_i(nblen(1),jpl) )
-               ALLOCATE( dta_bdy(ib_bdy)%h_i(nblen(1),jpl) )
-               ALLOCATE( dta_bdy(ib_bdy)%h_s(nblen(1),jpl) )
+               dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:)
+            ENDIF
+         ENDIF
+
+#if defined key_lim3
+         IF (cn_ice_lim(jbdy) /= 'none') THEN
+            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN
+               ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
+               ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
+               ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
             ELSE
-               IF ( ll_bdylim3 ) THEN ! case input is Xcat
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%a_i => bf(jfld)%fnow(:,1,:)
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%h_i => bf(jfld)%fnow(:,1,:)
-                  jfld = jfld + 1
-                  dta_bdy(ib_bdy)%h_s => bf(jfld)%fnow(:,1,:)
-               ELSE ! case input is 1cat
+               IF ( nice_cat == jpl ) THEN ! case input cat = jpl
+                  jfld = jfld + 1
+                  dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:)
+                  jfld = jfld + 1
+                  dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:)
+                  jfld = jfld + 1
+                  dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:)
+               ELSE                        ! case input cat = 1 OR (/=1 and /=jpl)
                   jfld_ai  = jfld + 1
                   jfld_hti = jfld + 2
                   jfld_hts = jfld + 3
                   jfld     = jfld + 3
-                  ALLOCATE( dta_bdy(ib_bdy)%a_i(nblen(1),jpl) )
-                  ALLOCATE( dta_bdy(ib_bdy)%h_i(nblen(1),jpl) )
-                  ALLOCATE( dta_bdy(ib_bdy)%h_s(nblen(1),jpl) )
-                  dta_bdy(ib_bdy)%a_i(:,:) = 0._wp
-                  dta_bdy(ib_bdy)%h_i(:,:) = 0._wp
-                  dta_bdy(ib_bdy)%h_s(:,:) = 0._wp
+                  ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) )
+                  ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) )
+                  ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) )
+                  dta_bdy(jbdy)%a_i(:,:) = 0._wp
+                  dta_bdy(jbdy)%h_i(:,:) = 0._wp
+                  dta_bdy(jbdy)%h_s(:,:) = 0._wp
                ENDIF
 
@@ -816 +823 @@
 #endif
          !
-      END DO ! ib_bdy 
+      END DO ! jbdy 
       !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_init')

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -48 +48 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-]
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice_tr     !: solar flux transmitted below the ice surface      [W/m2]
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s]
 
@@ -123 +125 @@
       ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
          &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
-         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   &
+         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl)  , alb_ice (jpi,jpj,jpl) ,   &
+         &      qml_ice(jpi,jpj,jpl)  , qcn_ice(jpi,jpj,jpl)   , qsr_ice_tr(jpi,jpj,jpl),  &
          &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   &
-         &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     &
          &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
          &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
@@ -139 +141 @@
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
                 STAT= ierr(2) )
-      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , &
-         &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , &
+      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , tn_ice (jpi,jpj,1)    , &
+         &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , &
          &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , &
          &                     STAT= ierr(3) )      
@@ -161 +163 @@
    PRIVATE
 
-   PUBLIC sbc_ice_alloc
+   PUBLIC   sbc_ice_alloc   !
 
    LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3 ice model
@@ -168 +170 @@
    REAL(wp)        , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-]
    INTEGER         , PUBLIC, PARAMETER ::   jpl = 1 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i
@@ -176 +178 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt, botmelt
    !
-   !! arrays relating to embedding ice in the ocean. These arrays need to be declared 
-   !! even if no ice model is required. In the no ice model or traditional levitating 
-   !! ice cases they contain only zeros
+   !! arrays related to embedding ice in the ocean. 
+   !! These arrays need to be declared even if no ice model is required. 
+   !! In the no ice model or traditional levitating ice cases they contain only zeros
    !! ---------------------
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_mass        !: mass of snow and ice at current  ice time step   [Kg/m2]

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

@@ -15 +15 @@
    !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
    !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
-   !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
+   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
    !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
+   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle) 
    !!----------------------------------------------------------------------
 
@@ -23 +24 @@
    !!   sbc_blk       : bulk formulation as ocean surface boundary condition
    !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
-   !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice
    !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
    !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
    !!   q_sat         : saturation humidity as a function of SLP and temperature
    !!   L_vap         : latent heat of vaporization of water as a function of temperature
+   !!             sea-ice case only : 
+   !!   blk_ice_tau   : provide the air-ice stress
+   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
+   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler)
+   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
+   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -42 +48 @@
    USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
    USE icethd_dh      ! for CALL ice_thd_snwblow
+   USE icethd_zdf,ONLY:   rn_cnd_s  ! for blk_ice_qcn
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
@@ -63 +70 @@
    PUBLIC   blk_ice_tau   ! routine called in icestp module
    PUBLIC   blk_ice_flx   ! routine called in icestp module
-#endif
+   PUBLIC   blk_ice_qcn   ! routine called in icestp module
+#endif 
 
 !!Lolo: should ultimately be moved in the module with all physical constants ?
@@ -111 +119 @@
    LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
    !
-!!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
@@ -139 +146 @@
       !!             ***  ROUTINE sbc_blk_alloc ***
       !!-------------------------------------------------------------------
-!!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc )
       ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
          &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
@@ -147 +153 @@
    END FUNCTION sbc_blk_alloc
 
+
    SUBROUTINE sbc_blk_init
       !!---------------------------------------------------------------------
@@ -154 +161 @@
       !!
       !! ** Method  : 
-      !!
-      !!      C A U T I O N : never mask the surface stress fields
-      !!                      the stress is assumed to be in the (i,j) mesh referential
-      !!
-      !! ** Action  :   
       !!
       !!----------------------------------------------------------------------
@@ -285 +287 @@
       !!
       !! ** Purpose :   provide at each time step the surface ocean fluxes
-      !!      (momentum, heat, freshwater and runoff)
+      !!              (momentum, heat, freshwater and runoff)
       !!
       !! ** Method  : (1) READ each fluxes in NetCDF files:
@@ -566 +568 @@
    END SUBROUTINE blk_oce
 
+
+
+   FUNCTION rho_air( ptak, pqa, pslp )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION rho_air  ***
+      !!
+      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
+      !!-------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
+      !!-------------------------------------------------------------------------------
+      !
+      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
+      !
+   END FUNCTION rho_air
+
+
+   FUNCTION cp_air( pqa )
+      !!-------------------------------------------------------------------------------
+      !!                           ***  FUNCTION cp_air  ***
+      !!
+      !! ** Purpose : Compute specific heat (Cp) of moist air
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!-------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
+      !!-------------------------------------------------------------------------------
+      !
+      Cp_air = Cp_dry + Cp_vap * pqa
+      !
+   END FUNCTION cp_air
+
+
+   FUNCTION q_sat( ptak, pslp )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION q_sat  ***
+      !!
+      !! ** Purpose : Specific humidity at saturation in [kg/kg]
+      !!              Based on accurate estimate of "e_sat"
+      !!              aka saturation water vapor (Goff, 1957)
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) ::   ze_sat, ztmp   ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            ztmp = rt0 / ptak(ji,jj)
+            !
+            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
+            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
+               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
+               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
+               !
+            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
+            !
+         END DO
+      END DO
+      !
+   END FUNCTION q_sat
+
+
+   FUNCTION gamma_moist( ptak, pqa )
+      !!----------------------------------------------------------------------------------
+      !!                           ***  FUNCTION gamma_moist  ***
+      !!
+      !! ** Purpose : Compute the moist adiabatic lapse-rate.
+      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
+      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
+      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
+      !
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) :: zrv, ziRT        ! local scalar
+      !!----------------------------------------------------------------------------------
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
+            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
+            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
+         END DO
+      END DO
+      !
+   END FUNCTION gamma_moist
+
+
+   FUNCTION L_vap( psst )
+      !!---------------------------------------------------------------------------------
+      !!                           ***  FUNCTION L_vap  ***
+      !!
+      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
+      !!
+      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
+      !!----------------------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
+      !!----------------------------------------------------------------------------------
+      !
+      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
+      !
+   END FUNCTION L_vap
+
 #if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                       ESIM sea-ice model
+   !!----------------------------------------------------------------------
+   !!   blk_ice_tau : provide the air-ice stress
+   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
+   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler)
+   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
+   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 
+   !!----------------------------------------------------------------------
 
    SUBROUTINE blk_ice_tau
@@ -688 +818 @@
 
 
-   SUBROUTINE blk_ice_flx( ptsu, palb )
+   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE blk_ice_flx  ***
       !!
-      !! ** Purpose :   provide the surface boundary condition over sea-ice
+      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
       !!
       !! ** Method  :   compute heat and freshwater exchanged
@@ -701 +831 @@
       !!---------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
       !!
@@ -707 +839 @@
       REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
       REAL(wp) ::   zztmp, z1_lsub           !   -      -
+      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice
+      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
@@ -717 +851 @@
       IF( ln_timing )   CALL timing_start('blk_ice_flx')
       !
-      !
-      ! local scalars
-      zcoef_dqlw   = 4.0 * 0.95 * Stef
-      zcoef_dqla   = -Ls * 11637800. * (-5897.8)
+      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars
+      zcoef_dqla = -Ls * 11637800. * (-5897.8)
       !
       zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
@@ -815 +947 @@
       END DO
 
-      !--------------------------------------------------------------------
-      ! FRACTIONs of net shortwave radiation which is not absorbed in the
-      ! thin surface layer and penetrates inside the ice cover
-      ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
-      !
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+      ! --- absorbed and transmitted shortwave radiation (W/m2) --- !
+      !
+      ! former coding was
+      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+
+      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- !
+      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value
+      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1
+
+      qsr_ice_tr(:,:,:) = 0._wp
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1
+               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
+               !
+               IF ( phs(ji,jj,jl) <= 0._wp ) THEN   ;    zfrqsr_tr_i  = zfr1 + zfac * zfr2
+               ELSE                                 ;    zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque
+               ENDIF
+               !
+               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation 
+               !
+            END DO
+         END DO
+      END DO
       !
       IF(ln_ctl) THEN
@@ -836 +989 @@
    END SUBROUTINE blk_ice_flx
    
-#endif
-
-   FUNCTION rho_air( ptak, pqa, pslp )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION rho_air  ***
-      !!
-      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
-      !!-------------------------------------------------------------------------------
-      !
-      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
-      !
-   END FUNCTION rho_air
-
-
-   FUNCTION cp_air( pqa )
-      !!-------------------------------------------------------------------------------
-      !!                           ***  FUNCTION cp_air  ***
-      !!
-      !! ** Purpose : Compute specific heat (Cp) of moist air
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!-------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
-      !!-------------------------------------------------------------------------------
-      !
-      Cp_air = Cp_dry + Cp_vap * pqa
-      !
-   END FUNCTION cp_air
-
-
-   FUNCTION q_sat( ptak, pslp )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION q_sat  ***
-      !!
-      !! ** Purpose : Specific humidity at saturation in [kg/kg]
-      !!              Based on accurate estimate of "e_sat"
-      !!              aka saturation water vapor (Goff, 1957)
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) ::   ze_sat, ztmp   ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            ztmp = rt0 / ptak(ji,jj)
-            !
-            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
-            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
-               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
-               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
+
+   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_qcn  ***
+      !!
+      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
+      !!                to force sea ice / snow thermodynamics
+      !!                in the case JULES coupler is emulated
+      !!                
+      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
+      !!                following the 0-layer Semtner (1976) approach
+      !!
+      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
+      !!              - qcn_ice : surface inner conduction flux (W/m2)
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   k_monocat   ! single-category option
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu        ! sea ice / snow surface temperature
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb         ! sea ice base temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs         ! snow thickness
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi         ! sea ice thickness
+      !
+      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
+      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
+      !
+      INTEGER  ::   ji, jj, jl           ! dummy loop indices
+      INTEGER  ::   iter                 ! local integer
+      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
+      REAL(wp) ::   zkeff_h, ztsu        !
+      REAL(wp) ::   zqc, zqnet           !
+      REAL(wp) ::   zhe, zqa0            !
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
+      !!---------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn')
+      
+      ! -------------------------------------!
+      !      I   Enhanced conduction factor  !
+      ! -------------------------------------!
+      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3)
+      ! Fichefet and Morales Maqueda, JGR 1997
+      !
+      zgfac(:,:,:) = 1._wp
+      
+      SELECT CASE ( k_monocat ) 
+      !
+      CASE ( 1 , 3 )
+         !
+         zfac    =   1._wp /  ( rn_cnd_s + rcdic )
+         zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon
+         zfac3   =   2._wp / zepsilon
+         !   
+         DO jl = 1, jpl                
+            DO jj = 1 , jpj
+               DO ji = 1, jpi
+                  !                                ! Effective thickness
+                  zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac
+                  !                                ! Enhanced conduction factor
+                  IF( zhe >=  zfac2 )  &
+                     zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) )
+               END DO
+            END DO
+         END DO
+         !      
+      END SELECT
+      
+      ! -------------------------------------------------------------!
+      !      II   Surface temperature and conduction flux            !
+      ! -------------------------------------------------------------!
+      !
+      zfac   =   rcdic * rn_cnd_s 
+      !                             ! ========================== !
+      DO jl = 1, jpl                !  Loop over ice categories  !
+         !                          ! ========================== !
+         DO jj = 1 , jpj
+            DO ji = 1, jpi
+               !                    ! Effective conductivity of the snow-ice system divided by thickness
+               zkeff_h =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) )
+               !                    ! Store initial surface temperature
+               ztsu    =   ptsu(ji,jj,jl)
+               !                    ! Net initial atmospheric heat flux
+               zqa0    =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)
                !
-            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa]
-            !
-         END DO
-      END DO
-      !
-   END FUNCTION q_sat
-
-
-   FUNCTION gamma_moist( ptak, pqa )
-      !!----------------------------------------------------------------------------------
-      !!                           ***  FUNCTION gamma_moist  ***
-      !!
-      !! ** Purpose : Compute the moist adiabatic lapse-rate.
-      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
-      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
-      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
-      !
-      INTEGER  ::   ji, jj         ! dummy loop indices
-      REAL(wp) :: zrv, ziRT        ! local scalar
-      !!----------------------------------------------------------------------------------
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
-            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
-            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) )
-         END DO
-      END DO
-      !
-   END FUNCTION gamma_moist
-
-
-   FUNCTION L_vap( psst )
-      !!---------------------------------------------------------------------------------
-      !!                           ***  FUNCTION L_vap  ***
-      !!
-      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
-      !!
-      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
-      !!----------------------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
-      !!----------------------------------------------------------------------------------
-      !
-      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
-      !
-   END FUNCTION L_vap
-
-#if defined key_lim3
+               DO iter = 1, nit     ! --- Iteration loop
+                   !                      ! Conduction heat flux through snow-ice system (>0 downwards)
+                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) )
+                   !                      ! Surface energy budget
+                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc
+                   !                      ! Temperature update
+                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )
+               END DO
+               !
+               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
+               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
+               !
+            END DO
+         END DO
+         !
+      END DO 
+      !      
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn')
+      !
+   END SUBROUTINE blk_ice_qcn
+   
 
    SUBROUTINE Cdn10_Lupkes2012( Cd )
@@ -961 +1098 @@
       !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
       !!
-      !! ** Purpose :    Recompute the ice-atm drag at 10m height to make
-      !!                 it dependent on edges at leads, melt ponds and flows.
+      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m 
+      !!                 to make it dependent on edges at leads, melt ponds and flows.
       !!                 After some approximations, this can be resumed to a dependency
       !!                 on ice concentration.
@@ -1012 +1149 @@
       !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
       !!
-      !! ** pUrpose :    1lternative turbulent transfert coefficients formulation
+      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
       !!                 between sea-ice and atmosphere with distinct momentum 
       !!                 and heat coefficients depending on sea-ice concentration 

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90

@@ -49 +49 @@
    PUBLIC ::   TURB_COARE   ! called by sbcblk.F90
 
-   !                                           !! COARE own values for given constants:
-   REAL(wp), PARAMETER ::   zi0     = 600.      ! scale height of the atmospheric boundary layer...1
-   REAL(wp), PARAMETER ::   Beta0   =   1.25    ! gustiness parameter
-   REAL(wp), PARAMETER ::   rctv0   =   0.608   ! constant to obtain virtual temperature...
+   !                                              !! COARE own values for given constants:
+   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer...
+   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter
+   REAL(wp), PARAMETER ::   rctv0   =   0.608_wp   ! constant to obtain virtual temperature...
 
    !!----------------------------------------------------------------------

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -971 +971 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      USE zdf_oce,  ONLY : ln_zdfswm
-
-      IMPLICIT NONE
-
-      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
-      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
-      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
+      USE zdf_oce,  ONLY :   ln_zdfswm
+      !
+      INTEGER, INTENT(in) ::   kt          ! ocean model time step index
+      INTEGER, INTENT(in) ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
+      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
       !!
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
@@ -1170 +1168 @@
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
          IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
-                                                                    .OR. srcv(jpr_hsig)%laction ) THEN
+            &                                                       .OR. srcv(jpr_hsig)%laction ) THEN
             CALL sbc_stokes()
          ENDIF
@@ -1525 +1523 @@
    
 
-   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist )
+   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
@@ -1576 +1574 @@
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
-      ! optional arguments, used only in 'mixed oce-ice' case
+      !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case
       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
       REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
       REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
-      !
-      INTEGER ::   jl         ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
+      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
+      !
+      INTEGER  ::   ji, jj, jl   ! dummy loop index
+      REAL(wp) ::   ztri         ! local scalar
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zevap_ice, zdevap_ice
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice    !!gm , zfrqsr_tr_i
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
-      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
-      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
-
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
-      ziceld(:,:) = 1. - picefr(:,:)
-      zcptn(:,:) = rcp * sst_m(:,:)
+      ziceld(:,:) = 1._wp - picefr(:,:)
+      zcptn (:,:) = rcp * sst_m(:,:)
       !
       !                                                      ! ========================= !
@@ -1622 +1618 @@
 #if defined key_lim3
       ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
-      zsnw(:,:) = 0._wp  ;  CALL ice_thd_snwblow( ziceld, zsnw )
+      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw )
       
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
@@ -1659 +1655 @@
          sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
          tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
-         DO jl=1,jpl
+         DO jl = 1, jpl
             evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
             devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
-         ENDDO
+         END DO
       ELSE
-         emp_tot(:,:) =         zemp_tot(:,:)
-         emp_ice(:,:) =         zemp_ice(:,:)
-         emp_oce(:,:) =         zemp_oce(:,:)     
-         sprecip(:,:) =         zsprecip(:,:)
-         tprecip(:,:) =         ztprecip(:,:)
-         DO jl=1,jpl
+         emp_tot(:,:) = zemp_tot(:,:)
+         emp_ice(:,:) = zemp_ice(:,:)
+         emp_oce(:,:) = zemp_oce(:,:)     
+         sprecip(:,:) = zsprecip(:,:)
+         tprecip(:,:) = ztprecip(:,:)
+         DO jl = 1, jpl
             evap_ice (:,:,jl) = zevap_ice (:,:)
             devap_ice(:,:,jl) = zdevap_ice(:,:)
-         ENDDO
+         END DO
       ENDIF
 
@@ -1691 +1687 @@
         fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)
       ENDIF
-
+      !
       IF( ln_mixcpl ) THEN
          emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
@@ -1703 +1699 @@
          tprecip(:,:) =                                  ztprecip(:,:)
       ENDIF
-
+      !
 #endif
+
       ! outputs
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
@@ -1730 +1727 @@
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
-            ENDDO
+            END DO
          ENDIF
       CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
@@ -1743 +1740 @@
          ELSE
             qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
-            ENDDO
+            END DO
          ENDIF
       CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
@@ -1766 +1763 @@
       !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
       zqns_oce = 0._wp
-      WHERE( ziceld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
+      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
 
       ! Heat content per unit mass of snow (J/kg)
@@ -1859 +1856 @@
          ELSE
             ! Set all category values equal for the moment
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
-            ENDDO
+            END DO
          ENDIF
          zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1868 +1865 @@
          zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
-            ENDDO
+            END DO
          ELSE
             qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
-            DO jl=1,jpl
+            DO jl = 1, jpl
                zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
-            ENDDO
+            END DO
          ENDIF
       CASE( 'mixed oce-ice' )
@@ -1890 +1887 @@
       IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
          zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
-         DO jl=1,jpl
+         DO jl = 1, jpl
             zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
-         ENDDO
+         END DO
       ENDIF
 
@@ -1908 +1905 @@
          qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
          qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
-         DO jl=1,jpl
+         DO jl = 1, jpl
             qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
-         ENDDO
+         END DO
       ELSE
          qsr_tot(:,:  ) = zqsr_tot(:,:  )
@@ -1946 +1943 @@
       END SELECT
 
-      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
-      ! Used for LIM3
-      ! Coupled case: since cloud cover is not received from atmosphere 
-      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
-      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-
-      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw )
-      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
-      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
+      !                                                      ! ========================= !
+      !                                                      !      Transmitted Qsr      !   [W/m2]
+      !                                                      ! ========================= !
+      SELECT CASE( nice_jules )
+      CASE( np_jules_OFF    )       !==  No Jules coupler  ==!
+         !
+!!gm         ! former coding was    
+!!gm         ! Coupled case: since cloud cover is not received from atmosphere 
+!!gm         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+!!gm         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+!!gm         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+!!gm         
+!!gm         ! to retrieve that coding, we needed to access h_i & h_s from here
+!!gm         ! we could even retrieve cloud fraction from the coupler
+!!gm         !
+!!gm         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter
+!!gm         !
+!!gm         DO jl = 1, jpl
+!!gm            DO jj = 1 , jpj
+!!gm               DO ji = 1, jpi
+!!gm                  !              !--- surface transmission parameter (Grenfell Maykut 77) --- !
+!!gm                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice 
+!!gm                  !
+!!gm                  !              ! --- influence of snow and thin ice --- !
+!!gm                  IF ( phs(ji,jj,jl) >= 0.0_wp )   zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque
+!!gm                  IF ( phi(ji,jj,jl) <= 0.1_wp )   zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation
+!!gm               END DO
+!!gm            END DO
+!!gm         END DO
+!!gm         !
+!!gm         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation 
+!!gm         !
+!!gm better coding of the above calculation:
+         !
+         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
+         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77)
+         !
+         qsr_ice_tr(:,:,:) = ztri * qsr_ice(:,:,:)
+         WHERE( phs(:,:,:) >= 0.0_wp )   qsr_ice_tr(:,:,:) = 0._wp            ! snow fully opaque
+         WHERE( phi(:,:,:) <= 0.1_wp )   qsr_ice_tr(:,:,:) = qsr_ice(:,:,:)   ! thin ice transmits all solar radiation
+!!gm end
+         !     
+      CASE( np_jules_ACTIVE )       !==  Jules coupler is active  ==!
+         !
+         !                    ! ===> here we must receive the qsr_ice_tr array from the coupler
+         !                           for now just assume zero (fully opaque ice)
+         qsr_ice_tr(:,:,:) = 0._wp
+         !
+      END SELECT
       !
       IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx')
@@ -2081 +2116 @@
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
          ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
-         DO jl=1,jpl
+         DO jl = 1, jpl
             ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
          END DO

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

@@ -710 +710 @@
       !                                         !   Prepare Coupling fields   !
       !                                         ! =========================== !
-
-! x and y comp of ice velocity
-
+      !
+      ! x and y comp of ice velocity
+      !
       CALL cice2nemo(uvel,u_ice,'F', -1. )
       CALL cice2nemo(vvel,v_ice,'F', -1. )
-
-! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
-
-! Snow and ice thicknesses (CO_2 and CO_3)
-
-      DO jl = 1,ncat
+      !
+      ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
+      !
+      ! Snow and ice thicknesses (CO_2 and CO_3)
+      !
+      DO jl = 1, ncat
          CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
          CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
-      ENDDO
+      END DO
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam')

branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

@@ -21 +21 @@
    !
    USE in_out_manager ! I/O manager
-   USE iom            ! I/O manager library
+   USE iom            ! I/O library
    USE fldread        ! read input field at current time step
    USE lbclnk         !
@@ -41 +41 @@
    REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  r1_hisf_tbl            !:1/thickness of tbl
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  ralpha                 !:proportion of bottom cell influenced by tbl 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  misfkt, misfkb         !:Level of ice shelf base
-
-   LOGICAL, PUBLIC ::   l_isfcpl = .false.       ! isf recieved from oasis
-
-   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K]
-   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s]
-   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3]
-   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C]
-   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg]
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+
+   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
+
+   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
+   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
+   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      !: volumic mass of ice shelf              [kg/m3]
+   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C]
+   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
 
 !: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
@@ -70 +70 @@
    
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , LOCEAN-IPSL (2015)
+   !! NEMO/OPA 4.0 , LOCEAN-IPSL (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -91 +91 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       !
-      INTEGER :: ji, jj, jk   ! loop index
-      INTEGER :: ikt, ikb     ! local integers
+      INTEGER ::   ji, jj, jk   ! loop index
+      INTEGER ::   ikt, ikb     ! local integers
       REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep) 
       REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d