New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
2021WP/KNL-01_Sibylle_RK3_stage1 (diff) – NEMO

Changes between Version 22 and Version 23 of 2021WP/KNL-01_Sibylle_RK3_stage1


Ignore:
Timestamp:
2022-01-03T18:32:34+01:00 (2 years ago)
Author:
techene
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • 2021WP/KNL-01_Sibylle_RK3_stage1

    v22 v23  
    434434 
    435435 
    436  
    437  
     436- finally we wrote additionnal comments and did some cleaning ! 
     437 
     438 
     439{{{#!diff 
     440Index: DEC_dbg_RK3/src/OCE/stprk3.F90 
     441=================================================================== 
     442--- DEC_dbg_RK3/src/OCE/stprk3.F90      (revision 15591) 
     443+++ DEC_dbg_RK3/src/OCE/stprk3.F90      (working copy) 
     444@@ -156,11 +156,13 @@ 
     445 !!st                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
     446                          CALL zdf_phy( kstp, Nbb, Nbb, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
     447 !!gm 
     448+!!st check with gm !!!!! 
     449+                         CALL eos( ts(:,:,:,:,Nbb), rhd, rhop, gdept_0(:,:,:) ) ! now in situ density for tramle slope computation 
     450       !  LATERAL  PHYSICS 
     451       ! 
     452       IF( l_ldfslp ) THEN                             ! slope of lateral mixing 
     453 !!gm gdep 
     454-                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density 
     455+!!st         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density 
     456  
     457          IF( ln_zps .AND. .NOT. ln_isfcav)                                    & 
     458             &            CALL zps_hde    ( kstp, Nbb, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     459Index: DEC_dbg_RK3/src/OCE/IOM/restart.F90 
     460=================================================================== 
     461--- DEC_dbg_RK3/src/OCE/IOM/restart.F90 (revision 15591) 
     462+++ DEC_dbg_RK3/src/OCE/IOM/restart.F90 (working copy) 
     463@@ -176,7 +176,7 @@ 
     464          CALL iom_rstput( kt, nitrst, numrow, 'vn'  , vv(:,:,:       ,Kmm) ) 
     465          CALL iom_rstput( kt, nitrst, numrow, 'tn'  , ts(:,:,:,jp_tem,Kmm) ) 
     466          CALL iom_rstput( kt, nitrst, numrow, 'sn'  , ts(:,:,:,jp_sal,Kmm) ) 
     467-         IF( .NOT.lk_SWE )   CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop ) 
     468+!!$         IF( .NOT.lk_SWE )   CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop ) 
     469 #endif 
     470       ENDIF 
     471  
     472@@ -290,6 +290,9 @@ 
     473       CALL iom_get( numror, jpdom_auto, 'sb'   , ts(:,:,:,jp_sal,Kbb) ) 
     474       CALL iom_get( numror, jpdom_auto, 'uu_b' , uu_b(:,:       ,Kbb), cd_type = 'U', psgn = -1._wp ) 
     475       CALL iom_get( numror, jpdom_auto, 'vv_b' , vv_b(:,:       ,Kbb), cd_type = 'V', psgn = -1._wp ) 
     476+      uu(:,:,:  ,Kmm) = uu(:,:,:  ,Kbb)         ! all now value set to before for initialisation (sbc_ssm_init) 
     477+      vv(:,:,:  ,Kmm) = vv(:,:,:  ,Kbb) 
     478+      ts(:,:,:,:,Kmm) = ts(:,:,:,:,Kbb) 
     479 #else 
     480       !                             !*  Read Kmm fields   (MLF only) 
     481       IF(lwp) WRITE(numout,*)    '           Kmm u, v and T-S fields read in the restart file' 
     482@@ -313,27 +316,17 @@ 
     483       ENDIF 
     484 #endif 
     485       ! 
     486-      IF( .NOT.lk_SWE ) THEN 
     487-         IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN 
     488-            CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop )   ! now    potential density 
     489-         ELSE 
     490-#if defined key_RK3 
     491-            CALL eos( ts, Kbb, rhop ) 
     492-#else 
     493-            CALL eos( ts, Kmm, rhop ) 
     494-#endif 
     495-!!#if defined key_qco 
     496-!!            ALLOCATE( zgdept(jpi,jpj,jpk) ) 
     497-!!            DO jk = 1, jpk 
     498-!!               zgdept(:,:,jk) = gdept(:,:,jk,Kmm) 
     499-!!            END DO 
     500-!!            CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, zgdept ) 
     501-!!            DEALLOCATE( zgdept ) 
     502-!!#else 
     503-!!            CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept(:,:,:,Kmm) ) 
     504-!!#endif 
     505-         ENDIF 
     506-      ENDIF 
     507+!!$      IF( .NOT.lk_SWE ) THEN 
     508+!!$         IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN 
     509+!!$            CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop )   ! now    potential density 
     510+!!$         ELSE 
     511+!!$#if defined key_RK3 
     512+!!$            CALL eos( ts, Kbb, rhop ) ! not rhop but rhd  
     513+!!$#else 
     514+!!$            CALL eos( ts, Kmm, rhop ) 
     515+!!$#endif 
     516+!!$         ENDIF 
     517+!!$      ENDIF 
     518       ! 
     519    END SUBROUTINE rst_read 
     520  
     521@@ -377,7 +370,7 @@ 
     522          CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb) ) 
     523          ! 
     524          !                                     !*  RK3: Set ssh at Kmm for AGRIF 
     525-         ssh(:,:,Kmm) = 0._wp 
     526+         ssh(:,:,Kmm) = ssh(:,:,Kbb) 
     527 #else 
     528          !                                     !*  MLF: Read ssh at Kmm 
     529          IF(lwp) WRITE(numout,*) 
     530Index: DEC_dbg_RK3/src/OCE/DYN/divhor.F90 
     531=================================================================== 
     532--- DEC_dbg_RK3/src/OCE/DYN/divhor.F90  (revision 15591) 
     533+++ DEC_dbg_RK3/src/OCE/DYN/divhor.F90  (working copy) 
     534@@ -83,7 +83,7 @@ 
     535       !  
     536       pe3divUh(:,:,:) = 0._wp    !!gm to be applied to the halos only 
     537       ! 
     538-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                                   
     539+      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 ) 
     540          hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * puu(ji  ,jj,jk)      & 
     541             &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk)      & 
     542             &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * pvv(ji,jj  ,jk)      & 
     543@@ -100,7 +100,7 @@ 
     544       ! 
     545       IF( ln_isf )   CALL isf_hdiv( kt, Kmm, hdiv )            !==  + ice-shelf mass exchange ==! 
     546       ! 
     547-      CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change) 
     548+      IF( nn_hls==1 )   CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change) 
     549       ! 
     550 !!gm Patch before suppression of hdiv from all modules that use it 
     551 !      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            !==  e3t * Horizontal divergence  ==! 
     552@@ -107,9 +107,9 @@ 
     553 !         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
     554 !      END_3D 
     555 !JC: over whole domain, and after lbclnk on hdiv to prevent from reproducibility issues 
     556-      DO jk=1, jpkm1 
     557-         pe3divUh(:,:,jk) = hdiv(:,:,jk) * e3t(:,:,jk,Kmm) 
     558-      END DO 
     559+      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 ) 
     560+         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
     561+      END_3D 
     562 !!gm end 
     563       ! 
     564       ! 
     565Index: DEC_dbg_RK3/src/OCE/DYN/sshwzv.F90 
     566=================================================================== 
     567--- DEC_dbg_RK3/src/OCE/DYN/sshwzv.F90  (revision 15591) 
     568+++ DEC_dbg_RK3/src/OCE/DYN/sshwzv.F90  (working copy) 
     569@@ -328,40 +328,41 @@ 
     570          DO jk = 1, jpkm1 
     571             ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 
     572             ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 
     573-            DO_2D( 0, 0, 0, 0 ) 
     574+            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )  
     575                zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 
     576             END_2D 
     577          END DO 
     578-         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
     579+         IF( nn_hls == 1)   CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
     580          !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
     581          !                             ! Same question holds for hdiv. Perhaps just for security 
     582-         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
     583-            ! computation of w 
     584-            pww(:,:,jk) = pww(:,:,jk+1) - (   ze3div(:,:,jk) + zhdiv(:,:,jk)   & 
     585-               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       & 
     586-               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk) 
     587-         END DO 
     588-         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 
     589+         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence 
     590+            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (   ze3div(ji,jj,jk) + zhdiv(ji,jj,jk)   & 
     591+                 &                            + r1_Dt * (  e3t(ji,jj,jk,Kaa)       & 
     592+                 &                                       - e3t(ji,jj,jk,Kbb) )   ) * tmask(ji,jj,jk) 
     593+         END_3D 
     594+         ! 
     595          DEALLOCATE( zhdiv )  
     596          !                                            !=================================! 
     597       ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==! 
     598          !                                            !=================================! 
     599-         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence 
     600-            pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk) 
     601-         END DO 
     602+         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence 
     603+            pww(ji,jj,jk) = pww(ji,jj,jk+1) - ze3div(ji,jj,jk)  
     604+         END_3D 
     605          !                                            !==========================================! 
     606       ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco') 
     607          !                                            !==========================================! 
     608-         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence 
     609-            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]] 
     610-            pww(:,:,jk) = pww(:,:,jk+1) - (  ze3div(:,:,jk)                            & 
     611-               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk) 
     612-         END DO 
     613+         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence 
     614+         !                                                              ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]] 
     615+            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  ze3div(ji,jj,jk)                            & 
     616+               &                               + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) )  ) * tmask(ji,jj,jk) 
     617+         END_3D 
     618       ENDIF 
     619  
     620       IF( ln_bdy ) THEN 
     621          DO jk = 1, jpkm1 
     622-            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 
     623+            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 
     624+               pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj) 
     625+            END_2D 
     626          END DO 
     627       ENDIF 
     628       ! 
     629Index: DEC_dbg_RK3/src/OCE/DYN/dynspg_ts.F90 
     630=================================================================== 
     631--- DEC_dbg_RK3/src/OCE/DYN/dynspg_ts.F90       (revision 15591) 
     632+++ DEC_dbg_RK3/src/OCE/DYN/dynspg_ts.F90       (working copy) 
     633@@ -802,7 +802,7 @@ 
     634          END_2D 
     635 # endif    
     636          ! 
     637-         CALL lbc_lnk_multi( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions 
     638+         CALL lbc_lnk( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions 
     639          ! 
     640       ENDIF 
     641       ! 
     642Index: DEC_dbg_RK3/src/OCE/stprk3_stg.F90 
     643=================================================================== 
     644--- DEC_dbg_RK3/src/OCE/stprk3_stg.F90  (revision 15591) 
     645+++ DEC_dbg_RK3/src/OCE/stprk3_stg.F90  (working copy) 
     646@@ -33,6 +33,8 @@ 
     647 # if defined key_agrif 
     648    USE agrif_oce_interp 
     649 # endif 
     650+!!st -dbg stp2d (r3u) 
     651+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     652     
     653    ! 
     654    USE prtctl         ! print control 
     655@@ -147,17 +149,27 @@ 
     656       !                     !==  ssh/h0 ratio at Kaa  ==!  
     657       ! 
     658       IF( .NOT.lk_linssh )   CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )   ! "after" ssh/h_0 ratio at t,u,v-column 
     659+!!      SELECT CASE( kstg ) 
     660+         ! 
     661+!!      CASE ( 3 ) 
     662+         !!st required at each stage for div hor loops 
     663+         CALL lbc_lnk( 'stprk3_stg', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, r3f(:,:), 'F', 1._wp ) 
     664+         ! 
     665+!!      END SELECT 
     666       ! 
     667       ! 
     668       !                     !==  advective velocity at Kmm  ==! 
     669       ! 
     670       !                                            !- horizontal components -!   (zaU,zaV)  
     671-      zub(:,:) = un_adv(:,:)*r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)    ! barotropic velocity correction 
     672-      zvb(:,:) = vn_adv(:,:)*r1_hv(:,:,Kmm) - vv_b(:,:,Kmm) 
     673-      DO jk = 1, jpkm1                                         ! horizontal advective velocity 
     674-         zaU(:,:,jk) = uu(:,:,jk,Kmm) + zub(:,:)*umask(:,:,jk) 
     675-         zaV(:,:,jk) = vv(:,:,jk,Kmm) + zvb(:,:)*vmask(:,:,jk) 
     676-      END DO 
     677+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
     678+         zub(ji,jj) = un_adv(ji,jj)*r1_hu(ji,jj,Kmm) - uu_b(ji,jj,Kmm)    ! barotropic velocity correction 
     679+         zvb(ji,jj) = vn_adv(ji,jj)*r1_hv(ji,jj,Kmm) - vv_b(ji,jj,Kmm) 
     680+      END_2D 
     681+      DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )           ! horizontal advective velocity 
     682+         zaU(ji,jj,jk) = uu(ji,jj,jk,Kmm) + zub(ji,jj)*umask(ji,jj,jk) 
     683+         zaV(ji,jj,jk) = vv(ji,jj,jk,Kmm) + zvb(ji,jj)*vmask(ji,jj,jk) 
     684+      END_3D 
     685+ 
     686       !                                            !- vertical components -!   ww 
     687                          CALL wzv  ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )     ! ww cross-level velocity 
     688  
     689@@ -223,9 +235,7 @@ 
     690          IF(.NOT. ln_trcadv_mus ) THEN 
     691             ! 
     692             DO jn = 1, jptra 
     693-               DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     694-               tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero 
     695-               END_3D 
     696+               tr(:,:,:,jn,Krhs) = 0._wp                              ! set tracer trends to zero !!st ::: required because of tra_adv new loops 
     697             END DO 
     698             !                                      !==  advection of passive tracers  ==! 
     699             rDt_trc = rDt 
     700@@ -253,9 +263,7 @@ 
     701          !                 !---------------! 
     702          ! 
     703          DO jn = 1, jptra 
     704-            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     705-            tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero 
     706-            END_3D 
     707+            tr(:,:,:,jn,Krhs) = 0._wp                              ! set tracer trends to zero !!st ::: required because of tra_adv new loops 
     708          END DO 
     709          !                                         !==  advection of passive tracers  ==! 
     710          rDt_trc = rDt 
     711@@ -393,12 +401,11 @@ 
     712       zvb(A2D(0)) = vv_b(A2D(0),Kaa) - SUM( e3v_0(A2D(0),:)*vv(A2D(0),:,Kaa), 3 ) * r1_hv_0(A2D(0)) 
     713       ! 
     714       DO jk = 1, jpkm1                                            ! corrected horizontal velocity 
     715-         uu(:,:,jk,Kaa) = uu(:,:,jk,Kaa) + zub(:,:)*umask(:,:,jk) 
     716-         vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk) 
     717+         uu(A2D(0),jk,Kaa) = uu(A2D(0),jk,Kaa) + zub(A2D(0))*umask(A2D(0),jk) 
     718+         vv(A2D(0),jk,Kaa) = vv(A2D(0),jk,Kaa) + zvb(A2D(0))*vmask(A2D(0),jk) 
     719       END DO 
     720-!!st pourquoi ne pas mettre A2D(0) ici ?  
     721-          
     722  
     723+ 
     724       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     725       ! Set boundary conditions 
     726       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     727@@ -408,9 +415,12 @@ 
     728             CALL Agrif_dyn( kstp, kstg ) 
     729 # endif 
     730       !                                              !* local domain boundaries  (T-point, unchanged sign) 
     731-      CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   & 
     732-         &                             , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. ) 
     733+      CALL lbc_lnk( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   & 
     734+         &                       , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. ) 
     735+      ! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy 
     736+      IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp_RK3_stg', avm_k, 'W', 1.0_wp ) 
     737       ! 
     738+      ! 
     739       !                                              !* BDY open boundaries 
     740       IF( ln_bdy )   THEN 
     741                                CALL bdy_tra( kstp, Kbb, ts,     Kaa ) 
     742Index: DEC_dbg_RK3/src/OCE/ZDF/zdfphy.F90 
     743=================================================================== 
     744--- DEC_dbg_RK3/src/OCE/ZDF/zdfphy.F90  (revision 15591) 
     745+++ DEC_dbg_RK3/src/OCE/ZDF/zdfphy.F90  (working copy) 
     746@@ -209,7 +209,7 @@ 
     747       ioptio = 0 
     748       IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF 
     749       IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init          ;   ENDIF 
     750-      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init( Kmm )   ;   ENDIF 
     751+      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init          ;   ENDIF 
     752       IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init          ;   ENDIF 
     753       IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init( Kmm )   ;   ENDIF 
     754       ! 
     755@@ -350,7 +350,7 @@ 
     756       IF( ln_zdfiwm )   CALL zdf_iwm( kt, Kmm, avm, avt, avs )   ! internal wave (de Lavergne et al 2017) 
     757  
     758       !                                         !* Lateral boundary conditions (sign unchanged) 
     759-      IF(nn_hls==1) THEN 
     760+      IF(nn_hls==1) THEN                                         ! if nn_hls==2 lbc_lnk done in stpxxx routines 
     761          IF( l_zdfsh2 ) THEN 
     762             CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp,   & 
     763                   &                 avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp ) 
     764Index: DEC_dbg_RK3/src/OCE/ZDF/zdftke.F90 
     765=================================================================== 
     766--- DEC_dbg_RK3/src/OCE/ZDF/zdftke.F90  (revision 15591) 
     767+++ DEC_dbg_RK3/src/OCE/ZDF/zdftke.F90  (working copy) 
     768@@ -698,7 +698,7 @@ 
     769    END SUBROUTINE tke_avn 
     770  
     771  
     772-   SUBROUTINE zdf_tke_init( Kmm ) 
     773+   SUBROUTINE zdf_tke_init 
     774       !!---------------------------------------------------------------------- 
     775       !!                  ***  ROUTINE zdf_tke_init  *** 
     776       !! 
     777@@ -714,7 +714,6 @@ 
     778       !!---------------------------------------------------------------------- 
     779       USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag 
     780       !! 
     781-      INTEGER, INTENT(in) ::   Kmm          ! time level index 
     782       INTEGER             ::   ji, jj, jk   ! dummy loop indices 
     783       INTEGER             ::   ios 
     784       !! 
     785Index: DEC_dbg_RK3/src/OCE/TRA/traqsr.F90 
     786=================================================================== 
     787--- DEC_dbg_RK3/src/OCE/TRA/traqsr.F90  (revision 15591) 
     788+++ DEC_dbg_RK3/src/OCE/TRA/traqsr.F90  (working copy) 
     789@@ -179,8 +179,8 @@ 
     790 #endif 
     791          END_3D 
     792          !                                                     !- sea-ice : store the 1st level attenuation coefficient 
     793-         WHERE( etot3(A2D(0),1) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - etot3(A2D(0),2) / etot3(A2D(0),1) 
     794-         ELSEWHERE                           ;   fraqsr_1lev(A2D(0)) = 1._wp 
     795+         WHERE( etot3(A2D(nn_hls),1) /= 0._wp )   ;   fraqsr_1lev(A2D(nn_hls)) = 1._wp - etot3(A2D(nn_hls),2) / etot3(A2D(nn_hls),1) 
     796+         ELSEWHERE                                ;   fraqsr_1lev(A2D(nn_hls)) = 1._wp 
     797          END WHERE 
     798          ! 
     799       END SELECT 
     800@@ -207,9 +207,8 @@ 
     801             &                             / e3t(ji,jj,jk,Kmm) 
     802       END_3D 
     803       ! 
     804-!!st7-2 
     805       ! sea-ice: store the 1st ocean level attenuation coefficient 
     806-      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
     807+      DO_2D( 0, 0, 0, 0 ) 
     808          IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 
     809          ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp 
     810          ENDIF 
     811@@ -217,7 +216,6 @@ 
     812       !                                       !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain 
     813       !                                       !                otherwise restartability and reproducibility are broken  
     814       CALL lbc_lnk( 'tra_qsr', fraqsr_1lev(:,:), 'T', 1._wp ) 
     815-!!st      CALL lbc_lnk( 'tra_qsr', qsr_hc(:,:,:), 'T', 1._wp ) 
     816       ! 
     817       IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution 
     818          ALLOCATE( zetot(A2D(nn_hls),jpk) ) 
     819Index: DEC_dbg_RK3/src/OCE/TRA/traadv_fct.F90 
     820=================================================================== 
     821--- DEC_dbg_RK3/src/OCE/TRA/traadv_fct.F90      (revision 15591) 
     822+++ DEC_dbg_RK3/src/OCE/TRA/traadv_fct.F90      (working copy) 
     823@@ -186,6 +186,7 @@ 
     824                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     825                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     826             !                               ! update and guess with monotonic sheme 
     827+!!st pt should be updated with a DO_3D( 0,0,0,0, 1, jpkm1 ) values outside the interior are wrong but unused. 
     828             pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   & 
     829                &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 
     830             zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 
     831Index: DEC_dbg_RK3/src/OCE/TRA/trasbc.F90 
     832=================================================================== 
     833--- DEC_dbg_RK3/src/OCE/TRA/trasbc.F90  (revision 15591) 
     834+++ DEC_dbg_RK3/src/OCE/TRA/trasbc.F90  (working copy) 
     835@@ -275,7 +275,7 @@ 
     836              
     837 !!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 
     838       IF( .NOT.ln_traqsr  .AND. kstg == 1) THEN     ! no solar radiation penetration 
     839-         DO_2D( 0, 0, 0, 0 ) 
     840+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
     841             qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)         ! total heat flux in qns 
     842             qsr(ji,jj) = 0._wp                           ! qsr set to zero 
     843          END_2D 
     844Index: DEC_dbg_RK3/src/OCE/TRD/trdglo.F90 
     845=================================================================== 
     846--- DEC_dbg_RK3/src/OCE/TRD/trdglo.F90  (revision 15591) 
     847+++ DEC_dbg_RK3/src/OCE/TRD/trdglo.F90  (working copy) 
     848@@ -201,7 +201,7 @@ 
     849          zkz  (:,:,:) = 0._wp 
     850          zkepe(:,:,:) = 0._wp 
     851     
     852-         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density 
     853+!!clem st dbg         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density 
     854  
     855          zcof = 0.5_wp / rho0             ! Density flux at w-point 
     856          zkz(:,:,1) = 0._wp 
     857Index: DEC_dbg_RK3/src/TOP/trcini.F90 
     858=================================================================== 
     859--- DEC_dbg_RK3/src/TOP/trcini.F90      (revision 15591) 
     860+++ DEC_dbg_RK3/src/TOP/trcini.F90      (working copy) 
     861@@ -246,13 +246,6 @@ 
     862       ! 
     863       tr(:,:,:,:,Kaa) = 0._wp 
     864       ! 
     865-      IF( ln_trcbc .AND. lltrcbc )  THEN  
     866-        CALL trc_bc_ini ( jptra, Kmm  )            ! set tracers Boundary Conditions 
     867-        CALL trc_bc     ( nit000, Kmm, tr, Kaa )   ! tracers: surface and lateral Boundary Conditions 
     868-      ENDIF 
     869-      ! 
     870-      IF( ln_trcais ) CALL trc_ais_ini   ! set tracers from Antarctic Ice Sheet 
     871-      ! 
     872       IF( ln_rsttr ) THEN              ! restart from a file 
     873         ! 
     874         CALL trc_rst_read( Kbb, Kmm ) 
     875Index: DEC_dbg_RK3/src/TOP/trcrst.F90 
     876=================================================================== 
     877--- DEC_dbg_RK3/src/TOP/trcrst.F90      (revision 15591) 
     878+++ DEC_dbg_RK3/src/TOP/trcrst.F90      (working copy) 
     879@@ -159,9 +159,6 @@ 
     880       ! prognostic variables  
     881       ! -------------------- 
     882 #if defined key_RK3 
     883-      DO jn = 1, jptra      ! MLF : After time step (before the swap) put in TRN 
     884-         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kaa) ) 
     885-      END DO 
     886       DO jn = 1, jptra      ! RK3 : After time step (before the swap) put in TRB 
     887          CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kaa) ) 
     888       END DO 
    438889=== Documentation updates 
    439890