Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcice_cice.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcice_cice.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcice_cice.F90	(revision 13710)
@@ -12,5 +12,9 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+# if ! defined key_qco
    USE domvvl
+# else
+   USE domqco
+# endif
    USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi
    USE in_out_manager  ! I/O manager
@@ -213,10 +217,10 @@
 ! T point to U point
 ! T point to V point
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
          fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
       END_2D
 
-      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.,  fr_iv , 'V', 1. )
+      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.0_wp,  fr_iv , 'V', 1.0_wp )
 
       ! set the snow+ice mass
@@ -233,9 +237,12 @@
 !!gm This should be put elsewhere....   (same remark for limsbc)
 !!gm especially here it is assumed zstar coordinate, but it can be ztilde....
+#if defined key_qco
+            IF( .NOT.ln_linssh )   CALL dom_qco_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column
+#else
             IF( .NOT.ln_linssh ) THEN
                !
                DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*r1_ht_0(:,:)*tmask(:,:,jk) )
+                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*r1_ht_0(:,:)*tmask(:,:,jk) )
                ENDDO
                e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb)
@@ -267,4 +274,5 @@
                END DO
             ENDIF
+#endif
          ENDIF
       ENDIF
@@ -304,5 +312,5 @@
 ! x comp of wind stress (CI_1)
 ! U point to F point
-         DO_2D_10_11
+         DO_2D( 1, 0, 1, 1 )
             ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
                                  + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
@@ -312,5 +320,5 @@
 ! y comp of wind stress (CI_2)
 ! V point to F point
-         DO_2D_11_10
+         DO_2D( 1, 1, 1, 0 )
             ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
                                  + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
@@ -327,5 +335,5 @@
             qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * rLsub
 ! End of temporary code
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                IF(fr_i(ji,jj).eq.0.0) THEN
                   DO jl=1,ncat
@@ -429,5 +437,5 @@
 ! x comp and y comp of surface ocean current
 ! U point to F point
-      DO_2D_10_11
+      DO_2D( 1, 0, 1, 1 )
          ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
       END_2D
@@ -435,5 +443,5 @@
 
 ! V point to F point
-      DO_2D_11_10
+      DO_2D( 1, 1, 1, 0 )
          ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
       END_2D
@@ -459,5 +467,5 @@
 ! x comp and y comp of sea surface slope (on F points)
 ! T point to F point
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  )) * r1_e1u(ji,jj  )    &
             &               + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1)  ) * fmask(ji,jj,1)
@@ -466,5 +474,5 @@
 
 ! T point to F point
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj)) * r1_e2v(ji  ,jj)    &
             &               + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj)  ) *  fmask(ji,jj,1)
@@ -495,8 +503,8 @@
       ss_iou(:,:)=0.0
 ! F point to U point
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
       END_2D
-      CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1. )
+      CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1.0_wp )
 
 ! y comp of ocean-ice stress 
@@ -505,8 +513,8 @@
 ! F point to V point
 
-      DO_2D_10_00
+      DO_2D( 1, 0, 0, 0 )
          ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
       END_2D
-      CALL lbc_lnk( 'sbcice_cice', ss_iov , 'V', -1. )
+      CALL lbc_lnk( 'sbcice_cice', ss_iov , 'V', -1.0_wp )
 
 ! x and y comps of surface stress
@@ -561,5 +569,5 @@
       fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
       
-      CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1., sfx , 'T', 1. )
+      CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1.0_wp, sfx , 'T', 1.0_wp )
 
 ! Solar penetrative radiation and non solar surface heat flux
@@ -587,7 +595,7 @@
 #endif
       qsr(:,:)=qsr(:,:)+ztmp1(:,:)
-      CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1. )
-
-      DO_2D_11_11
+      CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1.0_wp )
+
+      DO_2D( 1, 1, 1, 1 )
          nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
       END_2D
@@ -600,5 +608,5 @@
       qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
 
-      CALL lbc_lnk( 'sbcice_cice', qns , 'T', 1. )
+      CALL lbc_lnk( 'sbcice_cice', qns , 'T', 1.0_wp )
 
 ! Prepare for the following CICE time-step
@@ -613,10 +621,10 @@
 ! T point to U point
 ! T point to V point
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
          fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
       END_2D
 
-      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1., fr_iv , 'V', 1. )
+      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.0_wp, fr_iv , 'V', 1.0_wp )
 
       ! set the snow+ice mass
@@ -872,6 +880,6 @@
 !        pcg(:,:)=0.0
          DO jn=1,jpnij
-            DO jj=nldjt(jn),nlejt(jn)
-               DO ji=nldit(jn),nleit(jn)
+            DO jj=njs0all(jn),nje0all(jn)
+               DO ji=nis0all(jn),nie0all(jn)
                   png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
                ENDDO
@@ -973,5 +981,5 @@
 
       pn(:,:)=0.0
-      DO_2D_10_10
+      DO_2D( 1, 0, 1, 0 )
          pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
       END_2D
@@ -993,6 +1001,6 @@
          png(:,:,:)=0.0
          DO jn=1,jpnij
-            DO jj=nldjt(jn),nlejt(jn)
-               DO ji=nldit(jn),nleit(jn)
+            DO jj=njs0all(jn),nje0all(jn)
+               DO ji=nis0all(jn),nie0all(jn)
                   png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
                ENDDO
