src/ICE/icethd_pnd.F90

@@ -84,7 +84,7 @@ CONTAINS
       INTEGER ::   ji, jj, jl        ! loop indices
       !!-------------------------------------------------------------------
 
-      ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) )
+      ALLOCATE( diag_dvpn_mlt(A2D(0)) , diag_dvpn_lid(A2D(0)) , diag_dvpn_drn(A2D(0)) , diag_dvpn_rnf(A2D(0))  )
       ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) )
       !
       diag_dvpn_mlt (:,:) = 0._wp   ;   diag_dvpn_drn (:,:) = 0._wp
@@ -98,7 +98,7 @@ CONTAINS
       at_i(:,:) = SUM( a_i, dim=3 )
       !
       DO jl = 1, jpl
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
             IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN
                wfx_pnd  (ji,jj)    = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice
                a_ip     (ji,jj,jl) = 0._wp
@@ -115,7 +115,7 @@ CONTAINS
       !  Identify grid cells with ice
       !------------------------------
       npti = 0   ;   nptidx(:) = 0
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF( at_i(ji,jj) >= epsi10 ) THEN
             npti = npti + 1
             nptidx( npti ) = (jj - 1) * jpi + ji
@@ -137,6 +137,8 @@ CONTAINS
          END SELECT
       ENDIF
 
+      ! the following fields need to be updated in the halos (done in icethd): a_ip, v_ip, v_il, h_ip, h_il
+
       !------------------------------------
       !  Diagnostics
       !------------------------------------
@@ -529,7 +531,7 @@ CONTAINS
          zv_pnd  , &     ! volume of meltwater contributing to ponds
          zv_mlt          ! total amount of meltwater produced
 
-      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp_ini , &   !! total melt pond water available before redistribution and drainage
+      REAL(wp), DIMENSION(A2D(0)) ::   zvolp_ini , &   !! total melt pond water available before redistribution and drainage
                                         zvolp     , &   !! total melt pond water volume
                                         zvolp_res       !! remaining melt pond water available after drainage
 
@@ -589,7 +591,7 @@ CONTAINS
       zvolp(:,:) = 0._wp
 
       DO jl = 1, jpl
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
 
                IF ( a_i(ji,jj,jl) > epsi10 ) THEN
 
@@ -637,7 +639,7 @@ CONTAINS
 
       IF( ln_pnd_lids ) THEN
 
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
 
             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp_ini(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
 
@@ -764,7 +766,7 @@ CONTAINS
 
       DO jl = 1, jpl
 
-         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         DO_2D( 0, 0, 0, 0 )
 
 !              ! zap lids on small ponds
 !              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
@@ -826,7 +828,7 @@ CONTAINS
        !!
        !!------------------------------------------------------------------
 
-       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
+       REAL (wp), DIMENSION(A2D(0)), INTENT(INOUT) :: &
           zvolp,                                       &  ! total available pond water
           zdvolp                                          ! remaining meltwater after redistribution
 
@@ -865,10 +867,10 @@ CONTAINS
 
       INTEGER  ::   ji, jj, jk, jl                    ! loop indices
 
-       a_ip(:,:,:) = 0._wp
-       h_ip(:,:,:) = 0._wp
+       a_ip(A2D(0),:) = 0._wp
+       h_ip(A2D(0),:) = 0._wp
 
-       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+       DO_2D( 0, 0, 0, 0 )
 
              IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
 

src/ICE/iceupdate.F90

@@ -55,7 +55,7 @@ CONTAINS
       !!-------------------------------------------------------------------
       !!             ***  ROUTINE ice_update_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc )
+      ALLOCATE( utau_oce(A2D(0)), vtau_oce(A2D(0)), tmod_io(A2D(1)), STAT=ice_update_alloc )
       !
       CALL mpp_sum( 'iceupdate', ice_update_alloc )
       IF( ice_update_alloc /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_alloc: failed to allocate arrays' )
@@ -104,22 +104,29 @@ CONTAINS
       ! Net heat flux on top of the ice-ocean (W.m-2)
       !----------------------------------------------
       IF( ln_cndflx ) THEN   ! ice-atm interface = conduction (and melting) fluxes
-         qt_atm_oi(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:) + &
-            &             SUM( a_i_b * ( qcn_ice + qml_ice + qtr_ice_top ), dim=3 ) + qemp_ice(:,:)
+         DO_2D( 0, 0, 0, 0 )
+            qt_atm_oi(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * ( qns_oce(ji,jj) + qsr_oce(ji,jj) ) + qemp_oce(ji,jj) &
+               &               +     SUM( a_i_b(ji,jj,:) * ( qcn_ice(ji,jj,:) + qml_ice(ji,jj,:) + qtr_ice_top(ji,jj,:) ) ) &
+               &               + qemp_ice(ji,jj)
+         END_2D
       ELSE                   ! ice-atm interface = solar and non-solar fluxes
-         qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 
+         DO_2D( 0, 0, 0, 0 )
+            qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
+         END_2D
       ENDIF
 
       ! --- case we bypass ice thermodynamics --- !
       IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
-         qt_atm_oi  (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
-         qt_oce_ai  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
-         emp_ice    (:,:)   = 0._wp
-         qemp_ice   (:,:)   = 0._wp
-         qevap_ice  (:,:,:) = 0._wp
+         DO_2D( 0, 0, 0, 0 )
+            qt_atm_oi  (ji,jj)   = ( 1._wp - at_i_b(ji,jj) ) * ( qns_oce(ji,jj) + qsr_oce(ji,jj) ) + qemp_oce(ji,jj)
+            qt_oce_ai  (ji,jj)   = ( 1._wp - at_i_b(ji,jj) ) *   qns_oce(ji,jj)                    + qemp_oce(ji,jj)
+            emp_ice    (ji,jj)   = 0._wp
+            qemp_ice   (ji,jj)   = 0._wp
+            qevap_ice  (ji,jj,:) = 0._wp
+         END_2D
       ENDIF
 
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
 
          ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2)
          !---------------------------------------------------
@@ -183,20 +190,22 @@ CONTAINS
          snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step
          !                                               ! new mass per unit area
          snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) )
-         !                                               ! time evolution of snow+ice mass
-         snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice
 
       END_2D
+      CALL lbc_lnk( 'iceupdate', snwice_mass, 'T', 1.0_wp, snwice_mass_b, 'T', 1.0_wp ) ! needed for sshwzv and dynspg_ts (lbc on emp is done in sbcmod)
+      
+      ! time evolution of snow+ice mass
+      snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_Dt_ice
 
       ! Storing the transmitted variables
       !----------------------------------
       fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction
-      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature
+      tn_ice(:,:,:) = t_su(A2D(0),:)        ! Ice surface temperature
 
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice ) ! ice albedo
-
+      CALL ice_alb( ln_pnd_alb, t_su(A2D(0),:), h_i(A2D(0),:), h_s(A2D(0),:), a_ip_eff(:,:,:), h_ip(A2D(0),:), cloud_fra(:,:), & ! <<== in
+         &                                                                                                      alb_ice(:,:,:) ) ! ==>> out
       !
       IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
          CALL update_rst( 'WRITE', kt )
@@ -216,8 +225,8 @@ CONTAINS
       IF( iom_use('sfxopw'  ) )   CALL iom_put( 'sfxopw', sfx_opw * 1.e-03 )   ! salt flux from open water formation
       IF( iom_use('sfxdyn'  ) )   CALL iom_put( 'sfxdyn', sfx_dyn * 1.e-03 )   ! salt flux from ridging rafting
       IF( iom_use('sfxbri'  ) )   CALL iom_put( 'sfxbri', sfx_bri * 1.e-03 )   ! salt flux from brines
-      IF( iom_use('sfxres'  ) )   CALL iom_put( 'sfxres', sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes
       IF( iom_use('sfxsub'  ) )   CALL iom_put( 'sfxsub', sfx_sub * 1.e-03 )   ! salt flux from sublimation
+      IF( iom_use('sfxres'  ) )   CALL iom_put( 'sfxres', sfx_res(A2D(0)) * 1.e-03 )   ! salt flux from undiagnosed processes
 
       ! --- mass fluxes [kg/m2/s] --- !
       CALL iom_put( 'emp_oce', emp_oce )   ! emp over ocean (taking into account the snow blown away from the ice)
@@ -232,13 +241,13 @@ CONTAINS
       CALL iom_put( 'vfxsni'    , wfx_sni     )   ! mass flux from snow-ice formation
       CALL iom_put( 'vfxopw'    , wfx_opw     )   ! mass flux from growth in open water
       CALL iom_put( 'vfxdyn'    , wfx_dyn     )   ! mass flux from dynamics (ridging)
-      CALL iom_put( 'vfxres'    , wfx_res     )   ! mass flux from undiagnosed processes
       CALL iom_put( 'vfxpnd'    , wfx_pnd     )   ! mass flux from melt ponds
       CALL iom_put( 'vfxsub'    , wfx_ice_sub )   ! mass flux from ice sublimation (ice-atm.)
       CALL iom_put( 'vfxsub_err', wfx_err_sub )   ! "excess" of sublimation sent to ocean
+      CALL iom_put( 'vfxres'    , wfx_res(A2D(0))     )   ! mass flux from undiagnosed processes
 
       IF ( iom_use( 'vfxthin' ) ) THEN   ! mass flux from ice growth in open water + thin ice (<20cm) => comparable to observations
-         ALLOCATE( z2d(jpi,jpj) )
+         ALLOCATE( z2d(A2D(0)) )
          WHERE( hm_i(:,:) < 0.2 .AND. hm_i(:,:) > 0. ) ; z2d = wfx_bog
          ELSEWHERE                                     ; z2d = 0._wp
          END WHERE
@@ -264,8 +273,8 @@ CONTAINS
       IF( iom_use('qtr_ice_top') )   CALL iom_put( 'qtr_ice_top', SUM( qtr_ice_top * a_i_b, dim=3 )                          )   !     solar flux transmitted thru ice surface
       IF( iom_use('qt_oce'     ) )   CALL iom_put( 'qt_oce'     ,      ( qsr_oce + qns_oce ) * ( 1._wp - at_i_b ) + qemp_oce )
       IF( iom_use('qt_ice'     ) )   CALL iom_put( 'qt_ice'     , SUM( ( qns_ice + qsr_ice ) * a_i_b, dim=3 )     + qemp_ice )
-      IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( 'qt_oce_ai'  , qt_oce_ai * tmask(:,:,1)                                   )   ! total heat flux at the ocean   surface: interface oce-(ice+atm)
-      IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( 'qt_atm_oi'  , qt_atm_oi * tmask(:,:,1)                                   )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce)
+      IF( iom_use('qt_oce_ai'  ) )   CALL iom_put( 'qt_oce_ai'  , qt_oce_ai * smask0                                         )   ! total heat flux at the ocean   surface: interface oce-(ice+atm)
+      IF( iom_use('qt_atm_oi'  ) )   CALL iom_put( 'qt_atm_oi'  , qt_atm_oi * smask0                                         )   ! total heat flux at the oce-ice surface: interface atm-(ice+oce)
       IF( iom_use('qemp_oce'   ) )   CALL iom_put( 'qemp_oce'   , qemp_oce                                                   )   ! Downward Heat Flux from E-P over ocean
       IF( iom_use('qemp_ice'   ) )   CALL iom_put( 'qemp_ice'   , qemp_ice                                                   )   ! Downward Heat Flux from E-P over ice
 
@@ -282,9 +291,9 @@ CONTAINS
       ! heat fluxes associated with mass exchange (freeze/melt/precip...)
       CALL iom_put ('hfxthd'     , hfx_thd     )   !
       CALL iom_put ('hfxdyn'     , hfx_dyn     )   !
-      CALL iom_put ('hfxres'     , hfx_res     )   !
       CALL iom_put ('hfxsub'     , hfx_sub     )   !
       CALL iom_put ('hfxspr'     , hfx_spr     )   ! Heat flux from snow precip heat content
+      CALL iom_put ('hfxres'     , hfx_res(A2D(0))     )   !
 
       ! other heat fluxes
       IF( iom_use('hfxsensib'  ) )   CALL iom_put( 'hfxsensib'  ,      qsb_ice_bot * at_i_b         )   ! Sensible oceanic heat flux
@@ -314,7 +323,7 @@ CONTAINS
       !!
       !! ** Action  : * at each ice time step (every nn_fsbc time step):
       !!                - compute the modulus of ice-ocean relative velocity
-      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
+      !!                  (*rho*Cd) at T-point (C-grid)
       !!                      tmod_io = rhoco * | U_ice-U_oce |
       !!                - update the modulus of stress at ocean surface
       !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
@@ -325,19 +334,19 @@ CONTAINS
       !!
       !!    NB: - ice-ocean rotation angle no more allowed
       !!        - here we make an approximation: taum is only computed every ice time step
-      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids
-      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
+      !!          This avoids mutiple average to pass from U,V grids to T grids
+      !!          taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
       !!
-      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
+      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (T-pts) updated with ice-ocean fluxes
       !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
       !!---------------------------------------------------------------------
       INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
       !
       INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
-      REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
-      REAL(wp) ::   zflagi                          !   -      -
+      REAL(wp) ::   zutau_ice, zu_t, zmodt   ! local scalar
+      REAL(wp) ::   zvtau_ice, zv_t, zrhoco  !   -      -
+      REAL(wp) ::   zflagi                   !   -      -
       !!---------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('iceupdate')
 
@@ -350,46 +359,51 @@ CONTAINS
       zrhoco = rho0 * rn_cio
       !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
-         DO_2D( 0, 0, 0, 0 )                          !* update the modulus of stress at ocean surface (T-point)
+         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    !*  rhoco * |U_ice-U_oce| at T-point
+            !                                               ! 2*(U_ice-U_oce) at T-point
+            zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
+            zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
+            !                                               ! |U_ice-U_oce|^2
+            zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
+            !
+            tmod_io(ji,jj) = zrhoco * SQRT( zmodt )
+         END_2D
+         IF( nn_hls == 1 )   CALL lbc_lnk( 'iceupdate', tmod_io, 'T', 1._wp )
+         !
+         DO_2D( 0, 0, 0, 0 )                                !* save the air-ocean stresses at ice time-step
             !                                               ! 2*(U_ice-U_oce) at T-point
-            zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)
-            zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)
-            !                                              ! |U_ice-U_oce|^2
+            zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj) ! u_oce = ssu_m
+            zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) ! v_oce = ssv_m
+            !                                               ! |U_ice-U_oce|^2
             zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
             !                                               ! update the ocean stress modulus
             taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * zrhoco * zmodt
-            tmod_io(ji,jj) = zrhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
+            !
+            utau_oce(ji,jj) = utau(ji,jj)
+            vtau_oce(ji,jj) = vtau(ji,jj)
          END_2D
-         CALL lbc_lnk( 'iceupdate', taum, 'T', 1.0_wp, tmod_io, 'T', 1.0_wp )
-         !
-         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
-         vtau_oce(:,:) = vtau(:,:)
          !
       ENDIF
       !
       !                                      !==  every ocean time-step  ==!
       IF ( ln_drgice_imp ) THEN
          ! Save drag with right sign to update top drag in the ocean implicit friction
-         rCdU_ice(:,:) = -r1_rho0 * tmod_io(:,:) * at_i(:,:) * tmask(:,:,1)
+         DO_2D( 1, 1, 1, 1 )
+            rCdU_ice(ji,jj) = -r1_rho0 * tmod_io(ji,jj) * at_i(ji,jj) * tmask(ji,jj,1)
+         END_2D
          zflagi = 0._wp
       ELSE
          zflagi = 1._wp
       ENDIF
       !
       DO_2D( 0, 0, 0, 0 )                             !* update the stress WITHOUT an ice-ocean rotation angle
-         ! ice area at u and v-points
-         zat_u  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji+1,jj    ) * tmask(ji+1,jj  ,1) )  &
-            &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji+1,jj  ,1) )
-         zat_v  = ( at_i(ji,jj) * tmask(ji,jj,1) + at_i (ji  ,jj+1  ) * tmask(ji  ,jj+1,1) )  &
-            &     / MAX( 1.0_wp , tmask(ji,jj,1) + tmask(ji  ,jj+1,1) )
          !                                                   ! linearized quadratic drag formulation
-         zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
-         zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+         zutau_ice   = 0.5_wp * tmod_io(ji,jj) * ( u_ice(ji,jj) + u_ice(ji-1,jj) - pu_oce(ji,jj) - pu_oce(ji-1,jj) )
+         zvtau_ice   = 0.5_wp * tmod_io(ji,jj) * ( v_ice(ji,jj) + v_ice(ji,jj-1) - pv_oce(ji,jj) - pv_oce(ji,jj-1) )
          !                                                   ! stresses at the ocean surface
-         utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
-         vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
+         utau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * utau_oce(ji,jj) + at_i(ji,jj) * zutau_ice
+         vtau(ji,jj) = ( 1._wp - at_i(ji,jj) ) * vtau_oce(ji,jj) + at_i(ji,jj) * zvtau_ice
       END_2D
-      CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition
       !
       IF( ln_timing )   CALL timing_stop('iceupdate')
       !
@@ -446,12 +460,12 @@ CONTAINS
                CALL iom_get( numrir, jpdom_auto, 'snwice_mass_b', snwice_mass_b )
             ELSE                                     ! start from rest
                IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
-               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
+               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) + rhow * (vt_ip(:,:) + vt_il(:,:)) )
                snwice_mass_b(:,:) = snwice_mass(:,:)
             ENDIF
          ELSE                                   !* Start from rest
             IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
-            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) )
+            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) + rhow * (vt_ip(:,:) + vt_il(:,:)) )
             snwice_mass_b(:,:) = snwice_mass(:,:)
          ENDIF
          !

src/NST/agrif_all_update.F90

@@ -112,12 +112,10 @@ CONTAINS
       CALL dom_qco_zgr( Kbb_a, Kmm_a ) 
 #endif
 #if defined key_si3
-      CALL lbc_lnk( 'finalize_lbc_for_agrif',  a_i, 'T',1._wp,  v_i,'T',1._wp,                 &
-           &                                   v_s, 'T',1._wp, sv_i,'T',1._wp, oa_i,'T',1._wp, &
-           &                                   a_ip,'T',1._wp, v_ip,'T',1._wp, v_il,'T',1._wp )
-      CALL lbc_lnk( 'finalize_lbc_for_agrif', t_su,'T',1._wp )
-      CALL lbc_lnk( 'finalize_lbc_for_agrif',  e_s,'T',1._wp )
-      CALL lbc_lnk( 'finalize_lbc_for_agrif',  e_i,'T',1._wp )
+      CALL lbc_lnk( 'finalize_lbc_for_agrif', a_i, 'T',1._wp,  v_i,'T',1._wp,                 &
+           &                                  v_s, 'T',1._wp, sv_i,'T',1._wp, oa_i,'T',1._wp, &
+           &                                  a_ip,'T',1._wp, v_ip,'T',1._wp, v_il,'T',1._wp, t_su,'T',1._wp )
+      CALL lbc_lnk( 'finalize_lbc_for_agrif', e_i,'T',1._wp, e_s,'T',1._wp )
       CALL lbc_lnk( 'finalize_lbc_for_agrif', u_ice, 'U', -1._wp, v_ice, 'V', -1._wp )
 #endif
 #if defined key_top

src/NST/agrif_ice_interp.F90

@@ -59,12 +59,10 @@ CONTAINS
       Agrif_UseSpecialValue = .TRUE.
       CALL Agrif_init_variable(tra_iceini_id,procname=interp_tra_ice)
       !
-      CALL lbc_lnk( 'agrif_istate_ice',  a_i,'T',1._wp,  v_i,'T',1._wp, &
-               &         v_s,'T',1._wp, sv_i,'T',1._wp, oa_i,'T',1._wp, &
-               &        a_ip,'T',1._wp, v_ip,'T',1._wp, v_il,'T',1._wp )
-      CALL lbc_lnk( 'agrif_istate_ice', t_su,'T',1._wp )
-      CALL lbc_lnk( 'agrif_istate_ice',  e_s,'T',1._wp )
-      CALL lbc_lnk( 'agrif_istate_ice',  e_i,'T',1._wp )
+      CALL lbc_lnk( 'agrif_istate_ice', a_i,'T',1._wp,  v_i,'T',1._wp, &
+               &                        v_s,'T',1._wp, sv_i,'T',1._wp, oa_i,'T',1._wp, &
+               &                        a_ip,'T',1._wp, v_ip,'T',1._wp, v_il,'T',1._wp, t_su,'T',1._wp )
+      CALL lbc_lnk( 'agrif_istate_ice', e_i,'T',1._wp, e_s,'T',1._wp )
       !
       ! Set u_ice, v_ice:
       use_sign_north = .TRUE.

src/NST/agrif_oce_update.F90

@@ -1893,7 +1893,7 @@ CONTAINS
                      DO jk=k1,k2-1
                         IF  (ABS((ptab(ji,jj,jk)-e3u_0(ji,jj,jk))*umask(ji,jj,jk)).GE.1.e-6)  THEN 
                            kindic_agr = kindic_agr + 1 
-                           print *, 'erro u-pt', mig0(ji), mjg0(jj), jk, mbku(ji,jj), ikbot, ptab(ji,jj,jk), e3u_0(ji,jj,jk)
+                           PRINT *, 'erro u-pt', mig(ji,0), mjg(jj,0), jk, mbku(ji,jj), ikbot, ptab(ji,jj,jk), e3u_0(ji,jj,jk)
                         ENDIF
                      END DO
                   ENDIF
@@ -1933,7 +1933,7 @@ CONTAINS
                      DO jk=k1,k2-1
                         IF  (ABS((ptab(ji,jj,jk)-e3v_0(ji,jj,jk))*vmask(ji,jj,jk)).GE.1.e-6)  THEN 
                            kindic_agr = kindic_agr + 1 
-                           print *, 'erro v-pt', mig0(ji), mjg0(jj), mbkv(ji,jj), ptab(ji,jj,jk), e3v_0(ji,jj,jk)
+                           PRINT *, 'erro v-pt', mig(ji,0), mjg(jj,0), mbkv(ji,jj), ptab(ji,jj,jk), e3v_0(ji,jj,jk)
                         ENDIF
                      END DO
                   ENDIF

src/NST/agrif_user.F90

@@ -1095,8 +1095,8 @@
       !!----------------------------------------------------------------------
       !
       SELECT CASE( i )
-      CASE(1)        ;   indglob = mig(indloc)
-      CASE(2)        ;   indglob = mjg(indloc)
+      CASE(1)        ;   indglob = mig(indloc,nn_hls)
+      CASE(2)        ;   indglob = mjg(indloc,nn_hls)
       CASE DEFAULT   ;   indglob = indloc
       END SELECT
       !
@@ -1115,10 +1115,10 @@
       INTEGER, INTENT(out) :: jmin, jmax
       !!----------------------------------------------------------------------
       !
-      imin = mig( 1 )
-      jmin = mjg( 1 )
-      imax = mig(jpi)
-      jmax = mjg(jpj)
+      imin = mig( 1 ,nn_hls)
+      jmin = mjg( 1 ,nn_hls)
+      imax = mig(jpi,nn_hls)
+      jmax = mjg(jpj,nn_hls)
       ! 
    END SUBROUTINE Agrif_get_proc_info
 

src/OCE/BDY/bdyini.F90

@@ -491,10 +491,10 @@ CONTAINS
       ! Find lenght of boundaries and rim on local mpi domain
       !------------------------------------------------------
       !
-      iwe = mig(1)
-      ies = mig(jpi)
-      iso = mjg(1) 
-      ino = mjg(jpj) 
+      iwe = mig(  1,nn_hls)
+      ies = mig(jpi,nn_hls)
+      iso = mjg(  1,nn_hls) 
+      ino = mjg(jpj,nn_hls) 
       !
       DO ib_bdy = 1, nb_bdy
          DO igrd = 1, jpbgrd
@@ -554,8 +554,8 @@ CONTAINS
                      & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
                      !
                      icount = icount  + 1
-                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy) - mig(1) + 1   ! global to local indexes
-                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy) - mjg(1) + 1   ! global to local indexes
+                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy) - mig(1,nn_hls) + 1   ! global to local indexes
+                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy) - mjg(1,nn_hls) + 1   ! global to local indexes
                      idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
                      idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
                   ENDIF
@@ -1014,7 +1014,7 @@ CONTAINS
                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)
-                  IF(  mig0(ii) > 2 .AND. mig0(ii) < Ni0glo-2 .AND. mjg0(ij) > 2 .AND. mjg0(ij) < Nj0glo-2  ) THEN
+                  IF(  mig(ii,0) > 2 .AND. mig(ii,0) < Ni0glo-2 .AND. mjg(ij,0) > 2 .AND. mjg(ij,0) < Nj0glo-2  ) THEN
                      WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
                      CALL ctl_stop( ctmp1 )
                   END IF
@@ -1090,7 +1090,7 @@ CONTAINS
                ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
                IF( i_offset == 1 .and. zefl + zwfl == 2._wp ) THEN
                   icount = icount + 1
-                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
+                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii,nn_hls),mjg(ij,nn_hls)
                ELSE
                   ztmp(ii,ij) = -zwfl + zefl
                ENDIF
@@ -1130,7 +1130,7 @@ CONTAINS
                znfl = zmask(ii,ij+j_offset  )
                ! This error check only works if you are using the bdyXmask arrays (which are set to 0 on rims)
                IF( j_offset == 1 .and. znfl + zsfl == 2._wp ) THEN
-                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
+                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii,nn_hls),mjg(ij,nn_hls)
                   icount = icount + 1
                ELSE
                   ztmp(ii,ij) = -zsfl + znfl
@@ -1594,8 +1594,8 @@ CONTAINS
          ztestmask(1:2)=0.
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-              IF( mig0(ji) == jpiwob(ib) .AND. mjg0(jj) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+              IF( mig(ji,0) == jpiwob(ib) .AND. mjg(jj,0) == jpjwdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+              IF( mig(ji,0) == jpiwob(ib) .AND. mjg(jj,0) == jpjwft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
@@ -1630,8 +1630,8 @@ CONTAINS
          ztestmask(1:2)=0.
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-              IF( mig0(ji) == jpieob(ib)+1 .AND. mjg0(jj) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+              IF( mig(ji,0) == jpieob(ib)+1 .AND. mjg(jj,0) == jpjedt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+              IF( mig(ji,0) == jpieob(ib)+1 .AND. mjg(jj,0) == jpjeft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
@@ -1666,8 +1666,8 @@ CONTAINS
          ztestmask(1:2)=0.
          DO ji = 1, jpi
             DO jj = 1, jpj             
-              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-              IF( mjg0(jj) == jpjsob(ib) .AND. mig0(ji) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+              IF( mjg(jj,0) == jpjsob(ib) .AND. mig(ji,0) == jpisdt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+              IF( mjg(jj,0) == jpjsob(ib) .AND. mig(ji,0) == jpisft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
@@ -1688,8 +1688,8 @@ CONTAINS
          ztestmask(1:2)=0.
          DO ji = 1, jpi
             DO jj = 1, jpj             
-               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
-               IF( mjg0(jj) == jpjnob(ib)+1 .AND. mig0(ji) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
+               IF( mjg(jj,0) == jpjnob(ib)+1 .AND. mig(ji,0) == jpindt(ib) )   ztestmask(1) = tmask(ji,jj,1)
+               IF( mjg(jj,0) == jpjnob(ib)+1 .AND. mig(ji,0) == jpinft(ib) )   ztestmask(2) = tmask(ji,jj,1)  
             END DO
          END DO
          CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain

src/OCE/CRS/crsfld.F90

@@ -211,8 +211,8 @@ CONTAINS
       
       !  sbc fields  
       CALL crs_dom_ope( ssh(:,:,Kmm) , 'VOL', 'T', tmask, sshn_crs , p_e12=e1e2t, p_e3=ze3t           , psgn=1.0_wp )  
-      CALL crs_dom_ope( utau , 'SUM', 'U', umask, utau_crs , p_e12=e2u  , p_surf_crs=e2u_crs  , psgn=1.0_wp )
-      CALL crs_dom_ope( vtau , 'SUM', 'V', vmask, vtau_crs , p_e12=e1v  , p_surf_crs=e1v_crs  , psgn=1.0_wp )
+      CALL crs_dom_ope( utau , 'SUM', 'T', tmask, utau_crs , p_e12=e2t  , p_surf_crs=e2t_crs  , psgn=1.0_wp ) !clem tau: check psgn ??
+      CALL crs_dom_ope( vtau , 'SUM', 'T', tmask, vtau_crs , p_e12=e1t  , p_surf_crs=e1t_crs  , psgn=1.0_wp ) !
       CALL crs_dom_ope( wndm , 'SUM', 'T', tmask, wndm_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )
       CALL crs_dom_ope( rnf  , 'MAX', 'T', tmask, rnf_crs                                     , psgn=1.0_wp )
       CALL crs_dom_ope( qsr  , 'SUM', 'T', tmask, qsr_crs  , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0_wp )

src/OCE/DIA/diacfl.F90

@@ -51,19 +51,19 @@ CONTAINS
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
-      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
-      REAL(wp)                         ::   zCu_max, zCv_max, zCw_max        ! local scalars
-      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
-      LOGICAL , DIMENSION(jpi,jpj,jpk) ::   llmsk
+      INTEGER                         ::   ji, jj, jk                       ! dummy loop indices
+      REAL(wp)                        ::   zCu_max, zCv_max, zCw_max        ! local scalars
+      INTEGER , DIMENSION(3)          ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
+      LOGICAL , DIMENSION(A2D(0),jpk) ::   llmsk
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('dia_cfl')
       !
-      llmsk(     1:nn_hls,:,:) = .FALSE.   ! exclude halos from the checked region
-      llmsk(Nie0+1:   jpi,:,:) = .FALSE.
-      llmsk(:,     1:nn_hls,:) = .FALSE.
-      llmsk(:,Nje0+1:   jpj,:) = .FALSE.
+      !llmsk(     1:nn_hls,:,:) = .FALSE.   ! exclude halos from the checked region
+      !llmsk(Nie0+1:   jpi,:,:) = .FALSE.
+      !llmsk(:,     1:nn_hls,:) = .FALSE.
+      !llmsk(:,Nje0+1:   jpj,:) = .FALSE.
       !
       DO_3D( 0, 0, 0, 0, 1, jpk )      ! calculate Courant numbers
          zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * rDt / e1u  (ji,jj)      ! for i-direction