Changeset 12377 for NEMO/trunk/src/OCE/DIA

Timestamp:

2020-02-12T15:39:06+01:00 (4 years ago)

Author:

acc

Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:

NEMO/trunk

Files:

• . (modified) (1 prop)
• src/OCE/DIA/dia25h.F90 (modified) (9 diffs)
• src/OCE/DIA/diaar5.F90 (modified) (25 diffs)
• src/OCE/DIA/diacfl.F90 (modified) (4 diffs)
• src/OCE/DIA/diadct.F90 (modified) (17 diffs)
• src/OCE/DIA/diadetide.F90 (copied) (copied from NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diadetide.F90)
• src/OCE/DIA/diaharm.F90 (deleted)
• src/OCE/DIA/diahsb.F90 (modified) (23 diffs)
• src/OCE/DIA/diahth.F90 (modified) (15 diffs)
• src/OCE/DIA/diamlr.F90 (copied) (copied from NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diamlr.F90)
• src/OCE/DIA/diaptr.F90 (modified) (24 diffs)
• src/OCE/DIA/diawri.F90 (modified) (46 diffs)

NEMO/trunk

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev@HEAD      ext/AGRIF
+^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev@HEAD      ext/AGRIF
+^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL

NEMO/trunk/src/OCE/DIA/dia25h.F90

@@ -39 +39 @@
 CONTAINS
 
-   SUBROUTINE dia_25h_init 
+   SUBROUTINE dia_25h_init( Kbb )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE dia_25h_init  ***
@@ -47 +47 @@
       !! ** Method : Read namelist
       !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb       ! Time level index
+      !
       INTEGER ::   ios                 ! Local integer output status for namelist read
       INTEGER ::   ierror              ! Local integer for memory allocation
@@ -53 +55 @@
       !!----------------------------------------------------------------------
       !
-      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
       READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist' )
-      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
       READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist' )
@@ -95 +95 @@
       ! ------------------------- !
       cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem)
-      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal)
-      sshn_25h(:,:)   = sshb(:,:)
-      un_25h  (:,:,:) = ub  (:,:,:)
-      vn_25h  (:,:,:) = vb  (:,:,:)
+      tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kbb)
+      sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kbb)
+      sshn_25h(:,:)   = ssh(:,:,Kbb)
+      un_25h  (:,:,:) = uu  (:,:,:,Kbb)
+      vn_25h  (:,:,:) = vv  (:,:,:,Kbb)
       avt_25h (:,:,:) = avt (:,:,:)
       avm_25h (:,:,:) = avm (:,:,:)
@@ -116 +116 @@
 
 
-   SUBROUTINE dia_25h( kt )  
+   SUBROUTINE dia_25h( kt, Kmm )  
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_25h  ***
@@ -125 +125 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !!
       INTEGER ::   ji, jj, jk
@@ -150 +151 @@
       ! wn_25h could not be initialised in dia_25h_init, so we do it here instead
       IF( kt == nn_it000 ) THEN
-         wn_25h(:,:,:) = wn(:,:,:)
+         wn_25h(:,:,:) = ww(:,:,:)
       ENDIF
 
@@ -161 +162 @@
          ENDIF
 
-         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem)
-         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal)
-         sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:)
-         un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:)
-         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:)
-         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:)
+         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + ts (:,:,:,jp_tem,Kmm)
+         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + ts (:,:,:,jp_sal,Kmm)
+         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Kmm)
+         un_25h  (:,:,:)     = un_25h  (:,:,:) + uu  (:,:,:,Kmm)
+         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vv  (:,:,:,Kmm)
+         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + ww  (:,:,:)
          avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
          avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
@@ -245 +246 @@
          !
          ! After the write reset the values to cnt=1 and sum values equal current value 
-         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem)
-         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal)
-         sshn_25h(:,:)   = sshn(:,:)
-         un_25h  (:,:,:) = un  (:,:,:)
-         vn_25h  (:,:,:) = vn  (:,:,:)
-         wn_25h  (:,:,:) = wn  (:,:,:)
+         tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kmm)
+         sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kmm)
+         sshn_25h(:,:)   = ssh(:,:,Kmm)
+         un_25h  (:,:,:) = uu  (:,:,:,Kmm)
+         vn_25h  (:,:,:) = vv  (:,:,:,Kmm)
+         wn_25h  (:,:,:) = ww  (:,:,:)
          avt_25h (:,:,:) = avt (:,:,:)
          avm_25h (:,:,:) = avm (:,:,:)

NEMO/trunk/src/OCE/DIA/diaar5.F90

@@ -39 +39 @@
       
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -62 +62 @@
 
 
-   SUBROUTINE dia_ar5( kt )
+   SUBROUTINE dia_ar5( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ar5  ***
@@ -70 +70 @@
       !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
       !
       INTEGER  ::   ji, jj, jk, iks, ikb                      ! dummy loop arguments
@@ -89 +90 @@
          ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
-         zarea_ssh(:,:) = area(:,:) * sshn(:,:)
+         zarea_ssh(:,:) = area(:,:) * ssh(:,:,Kmm)
       ENDIF
       !
@@ -99 +100 @@
          zrhd(:,:,jpk) = 0._wp        ! ocean volume ; rhd is used as workspace
          DO jk = 1, jpkm1
-            zrhd(:,:,jk) = area(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+            zrhd(:,:,jk) = area(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
          CALL iom_put( 'volcello'  , zrhd(:,:,:)  )  ! WARNING not consistent with CMIP DR where volcello is at ca. 2000
-         CALL iom_put( 'masscello' , rau0 * e3t_n(:,:,:) * tmask(:,:,:) )  ! ocean mass
+         CALL iom_put( 'masscello' , rau0 * e3t(:,:,:,Kmm) * tmask(:,:,:) )  ! ocean mass
       ENDIF 
       !
       IF( iom_use( 'e3tb' ) )  THEN    ! bottom layer thickness
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikb = mbkt(ji,jj)
-               z2d(ji,jj) = e3t_n(ji,jj,ikb)
-            END DO
-         END DO
+         DO_2D_11_11
+            ikb = mbkt(ji,jj)
+            z2d(ji,jj) = e3t(ji,jj,ikb,Kmm)
+         END_2D
          CALL iom_put( 'e3tb', z2d )
       ENDIF 
@@ -122 +121 @@
          CALL iom_put( 'voltot', zvol               )
          CALL iom_put( 'sshtot', zvolssh / area_tot )
-         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
+         CALL iom_put( 'sshdyn', ssh(:,:,Kmm) - (zvolssh / area_tot) )
          !
       ENDIF
@@ -128 +127 @@
       IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN    
          !                     
-         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
+         ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh
          ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-         CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
+         CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity
          !
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
          DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)
          END DO
          IF( ln_linssh ) THEN
@@ -141 +140 @@
                   DO jj = 1, jpj
                      iks = mikt(ji,jj)
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,iks) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
             END IF
 !!gm
@@ -157 +156 @@
       
          !                                         ! steric sea surface height
-         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
+         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density
          zrhop(:,:,jpk) = 0._wp
          CALL iom_put( 'rhop', zrhop )
@@ -163 +162 @@
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
          DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)
          END DO
          IF( ln_linssh ) THEN
@@ -170 +169 @@
                   DO jj = 1,jpj
                      iks = mikt(ji,jj)
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,iks) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
             END IF
          END IF
@@ -183 +182 @@
          !                                         ! ocean bottom pressure
          zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
-         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
+         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + ssh(:,:,Kmm) + thick0(:,:) )
          CALL iom_put( 'botpres', zbotpres )
          !
@@ -191 +190 @@
           !                                         ! Mean density anomalie, temperature and salinity
           ztsn(:,:,:,:) = 0._wp                    ! ztsn(:,:,1,jp_tem/sal) is used here as 2D Workspace for temperature & salinity
-          DO jk = 1, jpkm1
-             DO jj = 1, jpj
-                DO ji = 1, jpi
-                   zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
-                   ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * tsn(ji,jj,jk,jp_tem)
-                   ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * tsn(ji,jj,jk,jp_sal)
-                ENDDO
-             ENDDO
-          ENDDO
+          DO_3D_11_11( 1, jpkm1 )
+             zztmp = area(ji,jj) * e3t(ji,jj,jk,Kmm)
+             ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zztmp * ts(ji,jj,jk,jp_tem,Kmm)
+             ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zztmp * ts(ji,jj,jk,jp_sal,Kmm)
+          END_3D
 
           IF( ln_linssh ) THEN
@@ -206 +201 @@
                   DO jj = 1, jpj
                      iks = mikt(ji,jj)
-                     ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * tsn(ji,jj,iks,jp_tem) 
-                     ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * tsn(ji,jj,iks,jp_sal) 
+                     ztsn(ji,jj,1,jp_tem) = ztsn(ji,jj,1,jp_tem) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_tem,Kmm) 
+                     ztsn(ji,jj,1,jp_sal) = ztsn(ji,jj,1,jp_sal) + zarea_ssh(ji,jj) * ts(ji,jj,iks,jp_sal,Kmm) 
                   END DO
                END DO
             ELSE
-               ztsn(:,:,1,jp_tem) = ztsn(:,:,1,jp_tem) + zarea_ssh(:,:) * tsn(:,:,1,jp_tem) 
-               ztsn(:,:,1,jp_sal) = ztsn(:,:,1,jp_sal) + zarea_ssh(:,:) * tsn(:,:,1,jp_sal) 
+               ztsn(:,:,1,jp_tem) = ztsn(:,:,1,jp_tem) + zarea_ssh(:,:) * ts(:,:,1,jp_tem,Kmm) 
+               ztsn(:,:,1,jp_sal) = ztsn(:,:,1,jp_sal) + zarea_ssh(:,:) * ts(:,:,1,jp_sal,Kmm) 
             END IF
          ENDIF
@@ -233 +228 @@
             ztpot(:,:,jpk) = 0._wp
             DO jk = 1, jpkm1
-               ztpot(:,:,jk) = eos_pt_from_ct( tsn(:,:,jk,jp_tem), tsn(:,:,jk,jp_sal) )
+               ztpot(:,:,jk) = eos_pt_from_ct( ts(:,:,jk,jp_tem,Kmm), ts(:,:,jk,jp_sal,Kmm) )
             END DO
             !
@@ -242 +237 @@
                z2d(:,:) = 0._wp
                DO jk = 1, jpkm1
-                 z2d(:,:) = z2d(:,:) + area(:,:) * e3t_n(:,:,jk) * ztpot(:,:,jk)
+                 z2d(:,:) = z2d(:,:) + area(:,:) * e3t(:,:,jk,Kmm) * ztpot(:,:,jk)
                END DO
                ztemp = glob_sum( 'diaar5', z2d(:,:)  ) 
@@ -255 +250 @@
              IF( iom_use( 'tosmint_pot') ) THEN
                z2d(:,:) = 0._wp
-               DO jk = 1, jpkm1
-                  DO jj = 1, jpj
-                     DO ji = 1, jpi   ! vector opt.
-                        z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  ztpot(ji,jj,jk)
-                     END DO
-                  END DO
-               END DO
+               DO_3D_11_11( 1, jpkm1 )
+                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t(ji,jj,jk,Kmm) *  ztpot(ji,jj,jk)
+               END_3D
                CALL iom_put( 'tosmint_pot', z2d ) 
             ENDIF
@@ -268 +259 @@
       ELSE       
          IF( iom_use('ssttot') ) THEN   ! Output sst in case we use EOS-80
-            zsst  = glob_sum( 'diaar5', area(:,:) * tsn(:,:,1,jp_tem) )
+            zsst  = glob_sum( 'diaar5', area(:,:) * ts(:,:,1,jp_tem,Kmm) )
             CALL iom_put('ssttot', zsst / area_tot )
          ENDIF
@@ -280 +271 @@
          zpe(:,:) = 0._wp
          IF( ln_zdfddm ) THEN
-            DO jk = 2, jpk
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     IF( rn2(ji,jj,jk) > 0._wp ) THEN
-                        zrw = ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk)
-                        !
-                        zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
-                        zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
-                        !
-                        zpe(ji, jj) = zpe(ji,jj)   &
-                           &        -  grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
-                           &                   - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
-                     ENDIF
-                  END DO
-               END DO
-             END DO
+            DO_3D_11_11( 2, jpk )
+               IF( rn2(ji,jj,jk) > 0._wp ) THEN
+                  zrw = ( gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm) ) / e3w(ji,jj,jk,Kmm)
+                  !
+                  zaw = rab_n(ji,jj,jk,jp_tem) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_tem)* zrw
+                  zbw = rab_n(ji,jj,jk,jp_sal) * (1. - zrw) + rab_n(ji,jj,jk-1,jp_sal)* zrw
+                  !
+                  zpe(ji, jj) = zpe(ji,jj)   &
+                     &        -  grav * (  avt(ji,jj,jk) * zaw * (ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )  &
+                     &                   - avs(ji,jj,jk) * zbw * (ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) )
+               ENDIF
+            END_3D
           ELSE
-            DO jk = 1, jpk
-               DO ji = 1, jpi
-                  DO jj = 1, jpj
-                     zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rau0 * e3w_n(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
+            DO_3D_11_11( 1, jpk )
+               zpe(ji,jj) = zpe(ji,jj) + avt(ji,jj,jk) * MIN(0._wp,rn2(ji,jj,jk)) * rau0 * e3w(ji,jj,jk,Kmm)
+            END_3D
          ENDIF
           CALL iom_put( 'tnpeo', zpe )
@@ -320 +303 @@
 
 
-   SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva ) 
+   SUBROUTINE dia_ar5_hst( ktra, cptr, puflx, pvflx ) 
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ar5_htr ***
@@ -329 +312 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: puflx  ! u-flux of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx  ! v-flux of advection/diffusion
       !
       INTEGER    ::  ji, jj, jk
@@ -336 +319 @@
 
     
-      z2d(:,:) = pua(:,:,1) 
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk) 
-            END DO
-         END DO
-       END DO
+      z2d(:,:) = puflx(:,:,1) 
+      DO_3D_00_00( 1, jpkm1 )
+         z2d(ji,jj) = z2d(ji,jj) + puflx(ji,jj,jk) 
+      END_3D
        CALL lbc_lnk( 'diaar5', z2d, 'U', -1. )
        IF( cptr == 'adv' ) THEN
@@ -354 +333 @@
        ENDIF
        !
-       z2d(:,:) = pva(:,:,1) 
-       DO jk = 1, jpkm1
-          DO jj = 2, jpjm1
-             DO ji = fs_2, fs_jpim1   ! vector opt.
-                z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk) 
-             END DO
-          END DO
-       END DO
+       z2d(:,:) = pvflx(:,:,1) 
+       DO_3D_00_00( 1, jpkm1 )
+          z2d(ji,jj) = z2d(ji,jj) + pvflx(ji,jj,jk) 
+       END_3D
        CALL lbc_lnk( 'diaar5', z2d, 'V', -1. )
        IF( cptr == 'adv' ) THEN
@@ -406 +381 @@
          zvol0 (:,:) = 0._wp
          thick0(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
-               DO ji = 1, jpi
-                  idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
-                  zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj)
-                  thick0(ji,jj) = thick0(ji,jj) +  idep    
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
+            zvol0 (ji,jj) = zvol0 (ji,jj) +  idep * area(ji,jj)
+            thick0(ji,jj) = thick0(ji,jj) +  idep    
+         END_3D
          vol0 = glob_sum( 'diaar5', zvol0 )
          DEALLOCATE( zvol0 )
@@ -428 +399 @@
             sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
             IF( ln_zps ) THEN               ! z-coord. partial steps
-               DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
-                  DO ji = 1, jpi
-                     ik = mbkt(ji,jj)
-                     IF( ik > 1 ) THEN
-                        zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                        sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
-                     ENDIF
-                  END DO
-               END DO
+               DO_2D_11_11
+                  ik = mbkt(ji,jj)
+                  IF( ik > 1 ) THEN
+                     zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                     sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
+                  ENDIF
+               END_2D
             ENDIF
             !

NEMO/trunk/src/OCE/DIA/diacfl.F90

@@ -33 +33 @@
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -41 +41 @@
 CONTAINS
 
-   SUBROUTINE dia_cfl ( kt )
+   SUBROUTINE dia_cfl ( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dia_cfl  ***
@@ -49 +49 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
@@ -64 +65 @@
       !
       !                
-      DO jk = 1, jpk       ! calculate Courant numbers
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
-               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
-               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction
-            END DO
-         END DO         
-      END DO
+      DO_3D_11_11( 1, jpk )
+         zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * z2dt / e1u  (ji,jj)      ! for i-direction
+         zCv_cfl(ji,jj,jk) = ABS( vv(ji,jj,jk,Kmm) ) * z2dt / e2v  (ji,jj)      ! for j-direction
+         zCw_cfl(ji,jj,jk) = ABS( ww(ji,jj,jk) ) * z2dt / e3w(ji,jj,jk,Kmm)   ! for k-direction
+      END_3D
       !
       ! write outputs

NEMO/trunk/src/OCE/DIA/diadct.F90

@@ -123 +123 @@
       !!---------------------------------------------------------------------
 
-     REWIND( numnam_ref )              ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections
      READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901)
 901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' )
 
-     REWIND( numnam_cfg )              ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections
      READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 )
 902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' )
@@ -175 +173 @@
  
  
-  SUBROUTINE dia_dct( kt )
+  SUBROUTINE dia_dct( kt, Kmm )
      !!---------------------------------------------------------------------
      !!               ***  ROUTINE diadct  ***  
@@ -192 +190 @@
      !!               Reinitialise all relevant arrays to zero 
      !!---------------------------------------------------------------------
-     INTEGER, INTENT(in) ::   kt
+     INTEGER, INTENT(in) ::   kt    ! ocean time step
+     INTEGER, INTENT(in) ::   Kmm   ! time level index
      !
      INTEGER ::   jsec              ! loop on sections
@@ -232 +231 @@
 
            !Compute transport through section  
-           CALL transport(secs(jsec),lldebug,jsec) 
+           CALL transport(Kmm,secs(jsec),lldebug,jsec) 
 
         ENDDO
@@ -246 +245 @@
            ! Sum over each class 
            DO jsec=1,nb_sec 
-              CALL dia_dct_sum(secs(jsec),jsec) 
+              CALL dia_dct_sum(Kmm,secs(jsec),jsec) 
            ENDDO 
 
@@ -558 +557 @@
 
 
-   SUBROUTINE transport(sec,ld_debug,jsec)
+   SUBROUTINE transport(Kmm,sec,ld_debug,jsec)
      !!-------------------------------------------------------------------------------------------
      !!                     ***  ROUTINE transport  ***
@@ -578 +577 @@
      !!
      !!-------------------------------------------------------------------------------------------
+     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index
      TYPE(SECTION),INTENT(INOUT) :: sec
      LOGICAL      ,INTENT(IN)    :: ld_debug
@@ -673 +673 @@
             SELECT CASE( sec%direction(jseg) )
                CASE(0,1) 
-                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
-                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
-                  zrhop = interp(k%I,k%J,jk,'V',rhop) 
-                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
-                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
+                  ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) ) 
+                  zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
+                  zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm)    ) * vmask(k%I,k%J,1) 
                CASE(2,3) 
-                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
-                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
-                  zrhop = interp(k%I,k%J,jk,'U',rhop) 
-                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
-                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+                  ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) ) 
+                  zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
+                  zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
                END SELECT 
                !
-               zdep= gdept_n(k%I,k%J,jk) 
+               zdep= gdept(k%I,k%J,jk,Kmm) 
   
                SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction 
                CASE(0,1)   
                   zumid=0._wp
-                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
+                  zvmid=isgnv*vv(k%I,k%J,jk,Kmm)*vmask(k%I,k%J,jk) 
                CASE(2,3) 
-                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
+                  zumid=isgnu*uu(k%I,k%J,jk,Kmm)*umask(k%I,k%J,jk) 
                   zvmid=0._wp
                END SELECT 
@@ -699 +699 @@
                !zTnorm=transport through one cell; 
                !velocity* cell's length * cell's thickness 
-               zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     & 
-                  &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk) 
+               zTnorm = zumid*e2u(k%I,k%J) * e3u(k%I,k%J,jk,Kmm)     & 
+                  &   + zvmid*e1v(k%I,k%J) * e3v(k%I,k%J,jk,Kmm) 
 
 !!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
@@ -765 +765 @@
 
 
-  SUBROUTINE dia_dct_sum(sec,jsec) 
+  SUBROUTINE dia_dct_sum(Kmm,sec,jsec) 
      !!------------------------------------------------------------- 
      !! Purpose: Average the transport over nn_dctwri time steps  
@@ -784 +784 @@
      !! 
      !!------------------------------------------------------------- 
+     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index
      TYPE(SECTION),INTENT(INOUT) :: sec 
      INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
@@ -845 +846 @@
               SELECT CASE( sec%direction(jseg) ) 
               CASE(0,1) 
-                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
-                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
-                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
-                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                 ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) ) 
+                 zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
+                 zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0) 
 
               CASE(2,3) 
-                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
-                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
-                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
-                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
-                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+                 ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) ) 
+                 zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
+                 zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                 zsshn =  0.5*( ssh(k%I,k%J,Kmm)    + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
               END SELECT 
  
-              zdep= gdept_n(k%I,k%J,jk) 
+              zdep= gdept(k%I,k%J,jk,Kmm) 
   
               !------------------------------- 
@@ -1101 +1102 @@
 
 
-   FUNCTION interp(ki, kj, kk, cd_point, ptab)
+   FUNCTION interp(Kmm, ki, kj, kk, cd_point, ptab)
   !!----------------------------------------------------------------------
   !!
@@ -1162 +1163 @@
   !!----------------------------------------------------------------------
   !*arguments
+  INTEGER, INTENT(IN)                          :: Kmm          ! time level index
   INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
   CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
@@ -1196 +1198 @@
   IF( ln_sco )THEN   ! s-coordinate case
 
-     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp 
-     zdep1 = gdept_n(ii1,ij1,kk) - zdepu
-     zdep2 = gdept_n(ii2,ij2,kk) - zdepu
+     zdepu = ( gdept(ii1,ij1,kk,Kmm) +  gdept(ii2,ij2,kk,Kmm) ) * 0.5_wp 
+     zdep1 = gdept(ii1,ij1,kk,Kmm) - zdepu
+     zdep2 = gdept(ii2,ij2,kk,Kmm) - zdepu
 
      ! weights
@@ -1210 +1212 @@
   ELSE       ! full step or partial step case 
 
-     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
-     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
-     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
+     ze3t  = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm) 
+     zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm)
+     zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm)
 
      IF(kk .NE. 1)THEN
@@ -1245 +1247 @@
       IMPLICIT NONE
    END SUBROUTINE dia_dct_init
-   SUBROUTINE dia_dct( kt )
+
+   SUBROUTINE dia_dct( kt, Kmm )         ! Dummy routine
       IMPLICIT NONE
-      INTEGER, INTENT(in) ::   kt
+      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
+      INTEGER, INTENT( in ) :: Kmm  ! ocean time level index
+      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
    END SUBROUTINE dia_dct
    !

NEMO/trunk/src/OCE/DIA/diahsb.F90

@@ -17 +17 @@
    USE phycst         ! physical constants
    USE sbc_oce        ! surface thermohaline fluxes
+   USE isf_oce        ! ice shelf fluxes
    USE sbcrnf         ! river runoff
-   USE sbcisf         ! ice shelves
    USE domvvl         ! vertical scale factors
    USE traqsr         ! penetrative solar radiation
@@ -48 +48 @@
    REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   tmask_ini
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -58 +57 @@
 CONTAINS
 
-   SUBROUTINE dia_hsb( kt )
+   SUBROUTINE dia_hsb( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE dia_hsb  ***
@@ -69 +68 @@
       !!
       !!---------------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER    ::   ji, jj, jk                  ! dummy loop indice
@@ -86 +86 @@
       IF( ln_timing )   CALL timing_start('dia_hsb')      
       !
-      tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ;
-      tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ;
+      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
+      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
       ! ------------------------- !
       ! 1 - Trends due to forcing !
       ! ------------------------- !
-      z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) )   ! volume fluxes
+      z_frc_trd_v = r1_rau0 * glob_sum( 'diahsb', - ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * surf(:,:) )   ! volume fluxes
       z_frc_trd_t =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes
       z_frc_trd_s =           glob_sum( 'diahsb', sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes
@@ -98 +98 @@
       IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( 'diahsb', rnf_tsc(:,:,jp_sal) * surf(:,:) )
       !                    ! Add ice shelf heat & salt input
-      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( 'diahsb', risf_tsc(:,:,jp_tem) * surf(:,:) )
+      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t &
+         &                          + glob_sum( 'diahsb', ( risf_cav_tsc(:,:,jp_tem) + risf_par_tsc(:,:,jp_tem) ) * surf(:,:) )
       !                    ! Add penetrative solar radiation
       IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( 'diahsb', qsr     (:,:) * surf(:,:) )
@@ -108 +109 @@
             DO ji=1,jpi
                DO jj=1,jpj
-                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem)
-                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal)
+                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
                END DO
             END DO
          ELSE
-            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem)
-            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal)
+            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
          END IF
          z_wn_trd_t = - glob_sum( 'diahsb', z2d0 ) 
@@ -135 +136 @@
 
       !                    ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*sshn(:,:) - surf_ini(:,:)*ssh_ini(:,:) )
+      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
 
       !                    ! heat & salt content variation (associated with ssh)
@@ -142 +143 @@
             DO ji = 1, jpi
                DO jj = 1, jpj
-                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
-                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) ) 
+                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) ) 
                END DO
             END DO
          ELSE                          ! no under ice-shelf seas
-            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 
-            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
+            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) ) 
+            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) ) 
          END IF
          z_ssh_hc = glob_sum_full( 'diahsb', z2d0 ) 
@@ -155 +156 @@
       !
       DO jk = 1, jpkm1           ! volume variation (calculated with scale factors)
-         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk) ) * tmask(:,:,jk)
+         zwrk(:,:,jk) = surf(:,:)*e3t(:,:,jk,Kmm)*tmask(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk)*tmask_ini(:,:,jk)
       END DO
-      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+      zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )     ! glob_sum_full needed as tmask and tmask_ini could be different
       DO jk = 1, jpkm1           ! heat content variation
-         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_tem) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
+         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) - surf_ini(:,:)*hc_loc_ini(:,:,jk) )
       END DO
       zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
       DO jk = 1, jpkm1           ! salt content variation
-         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_sal) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
+         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) - surf_ini(:,:)*sc_loc_ini(:,:,jk) )
       END DO
       zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
@@ -185 +186 @@
       ! ----------------------- !
       DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors)
-         zwrk(:,:,jk) = surf(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+         zwrk(:,:,jk) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
-      zvol_tot = glob_sum_full( 'diahsb', zwrk(:,:,:) )
+      zvol_tot = glob_sum( 'diahsb', zwrk(:,:,:) )
 
 !!gm to be added ?
 !      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution
-!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
+!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * ssh(:,:,Kmm) )
 !      ENDIF
 !!gm end
@@ -233 +234 @@
       ENDIF
       !
-      IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
+      IF( lrst_oce )   CALL dia_hsb_rst( kt, Kmm, 'WRITE' )
       !
       IF( ln_timing )   CALL timing_stop('dia_hsb')
@@ -240 +241 @@
 
 
-   SUBROUTINE dia_hsb_rst( kt, cdrw )
+   SUBROUTINE dia_hsb_rst( kt, Kmm, cdrw )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE dia_hsb_rst  ***
@@ -249 +250 @@
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index
       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !
@@ -270 +272 @@
             CALL iom_get( numror, jpdom_autoglo, 'ssh_ini'   , ssh_ini   , ldxios = lrxios )
             CALL iom_get( numror, jpdom_autoglo, 'e3t_ini'   , e3t_ini   , ldxios = lrxios )
+            CALL iom_get( numror, jpdom_autoglo, 'tmask_ini' , tmask_ini , ldxios = lrxios )
             CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini, ldxios = lrxios )
             CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini, ldxios = lrxios )
@@ -281 +284 @@
             IF(lwp) WRITE(numout,*)
             surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
-            ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
+            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh
             DO jk = 1, jpk
               ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
-               e3t_ini   (:,:,jk) = e3t_n(:,:,jk)                      * tmask(:,:,jk)  ! initial vertical scale factors
-               hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content
-               sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content
+               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors
+               tmask_ini (:,:,jk) = tmask(:,:,jk)                                       ! initial mask
+               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content
+               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content
             END DO
             frc_v = 0._wp                                           ! volume       trend due to forcing
@@ -295 +299 @@
                   DO ji = 1, jpi
                      DO jj = 1, jpj
-                        ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
-                        ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh
                      END DO
                    END DO
                 ELSE
-                  ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
-                  ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm)   ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm)   ! initial salt content in ssh
                END IF
                frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
@@ -325 +329 @@
          CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini'   , ssh_ini   , ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini'   , e3t_ini   , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'tmask_ini' , tmask_ini , ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini, ldxios = lwxios )
          CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini, ldxios = lwxios )
@@ -338 +343 @@
 
 
-   SUBROUTINE dia_hsb_init
+   SUBROUTINE dia_hsb_init( Kmm )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE dia_hsb  ***
@@ -350 +355 @@
       !!             - Compute coefficients for conversion
       !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! time level index
+      !
       INTEGER ::   ierror, ios   ! local integer
       !!
@@ -360 +367 @@
          WRITE(numout,*) '~~~~~~~~~~~~ '
       ENDIF
-      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
       READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' )
-      REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
       READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
@@ -396 +401 @@
       ! ------------------- !
       ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
-         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  )
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), tmask_ini(jpi,jpj,jpk),STAT=ierror  )
       IF( ierror > 0 ) THEN
          CALL ctl_stop( 'dia_hsb_init: unable to allocate hc_loc_ini' )   ;   RETURN
@@ -417 +422 @@
       ! 4 - initial conservation variables !
       ! ---------------------------------- !
-      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
+      CALL dia_hsb_rst( nit000, Kmm, 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE dia_hsb_init

NEMO/trunk/src/OCE/DIA/diahth.F90

@@ -40 +40 @@
 
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -61 +63 @@
 
 
-   SUBROUTINE dia_hth( kt )
+   SUBROUTINE dia_hth( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_hth  ***
@@ -82 +84 @@
       !!-------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
       !!
       INTEGER                      ::   ji, jj, jk            ! dummy loop arguments
@@ -126 +129 @@
             zdepinv(:,:) = 0._wp  
             zmaxdzT(:,:) = 0._wp  
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
-                  hth     (ji,jj) = zztmp
-                  zabs2   (ji,jj) = zztmp
-                  ztm2    (ji,jj) = zztmp
-                  zrho10_3(ji,jj) = zztmp
-                  zpycn   (ji,jj) = zztmp
-                 END DO
-            END DO
+            DO_2D_11_11
+               zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
+               hth     (ji,jj) = zztmp
+               zabs2   (ji,jj) = zztmp
+               ztm2    (ji,jj) = zztmp
+               zrho10_3(ji,jj) = zztmp
+               zpycn   (ji,jj) = zztmp
+            END_2D
             IF( nla10 > 1 ) THEN 
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
-                     zrho0_3(ji,jj) = zztmp
-                     zrho0_1(ji,jj) = zztmp
-                  END DO
-               END DO
+               DO_2D_11_11
+                  zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
+                  zrho0_3(ji,jj) = zztmp
+                  zrho0_1(ji,jj) = zztmp
+               END_2D
             ENDIF
       
             ! Preliminary computation
             ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( tmask(ji,jj,nla10) == 1. ) THEN
-                     zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)  &
-                        &           - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
-                        &           - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
-                     zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)  &
-                        &           - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
-                     zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
-                     zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
-                     zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
-                     zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
-                  ELSE
-                     zdelr(ji,jj) = 0._wp
-                  ENDIF
-               END DO
-            END DO
+            DO_2D_11_11
+               IF( tmask(ji,jj,nla10) == 1. ) THEN
+                  zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)  &
+                     &           - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   &
+                     &           - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
+                  zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)  &
+                     &           - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
+                  zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm)
+                  zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
+                  zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
+                  zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
+               ELSE
+                  zdelr(ji,jj) = 0._wp
+               ENDIF
+            END_2D
 
             ! ------------------------------------------------------------- !
@@ -172 +169 @@
             ! MLD: rho = rho(1) + zrho1                                     !
             ! ------------------------------------------------------------- !
-            DO jk = jpkm1, 2, -1   ! loop from bottom to 2
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     zzdep = gdepw_n(ji,jj,jk)
-                     zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) &
-                            & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
-                     zzdep = zzdep * tmask(ji,jj,1)
-
-                     IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
-                         zmaxdzT(ji,jj) = zztmp   
-                         hth    (ji,jj) = zzdep                ! max and depth of dT/dz
-                     ENDIF
-               
-                     IF( nla10 > 1 ) THEN 
-                        zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1)
-                        IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03
-                        IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01
-                     ENDIF
-                  END DO
-               END DO
-            END DO
+            DO_3DS_11_11( jpkm1, 2, -1 )
+               !
+               zzdep = gdepw(ji,jj,jk,Kmm)
+               zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) &
+                      & / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
+               zzdep = zzdep * tmask(ji,jj,1)
+
+               IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
+                   zmaxdzT(ji,jj) = zztmp   
+                   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
+               ENDIF
+         
+               IF( nla10 > 1 ) THEN 
+                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                       ! delta rho(1)
+                  IF( zztmp > zrho3 )   zrho0_3(ji,jj) = zzdep                ! > 0.03
+                  IF( zztmp > zrho1 )   zrho0_1(ji,jj) = zzdep                ! > 0.01
+               ENDIF
+            END_3D
          
             CALL iom_put( 'mlddzt', hth )            ! depth of the thermocline
@@ -213 +206 @@
             ! depth of temperature inversion                                !
             ! ------------------------------------------------------------- !
-            DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
-                     !
-                     zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
-                     IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
-                     IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
-                     zztmp = -zztmp                                          ! delta T(10m)
-                     IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
-                        ztinv(ji,jj) = zztmp   
-                        zdepinv (ji,jj) = zzdep   ! max value and depth
-                     ENDIF
-
-                     zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
-                     IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
-                     IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
-                     !
-                  END DO
-               END DO
-            END DO
+            DO_3DS_11_11( jpkm1, nlb10, -1 )
+               !
+               zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
+               !
+               zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm)  ! - delta T(10m)
+               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
+               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
+               zztmp = -zztmp                                          ! delta T(10m)
+               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
+                  ztinv(ji,jj) = zztmp   
+                  zdepinv (ji,jj) = zzdep   ! max value and depth
+               ENDIF
+
+               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
+               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
+               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
+               !
+            END_3D
 
             CALL iom_put( 'mld_dt02', zabs2    )   ! MLD abs(delta t) - 0.2
@@ -250 +239 @@
          IF( iom_use ('20d') ) THEN  ! depth of the 20 isotherm
             ztem2 = 20.
-            CALL dia_hth_dep( ztem2, hd20 )  
+            CALL dia_hth_dep( Kmm, ztem2, hd20 )  
             CALL iom_put( '20d', hd20 )    
          ENDIF
@@ -256 +245 @@
          IF( iom_use ('26d') ) THEN  ! depth of the 26 isotherm
             ztem2 = 26.
-            CALL dia_hth_dep( ztem2, hd26 )  
+            CALL dia_hth_dep( Kmm, ztem2, hd26 )  
             CALL iom_put( '26d', hd26 )    
          ENDIF
@@ -262 +251 @@
          IF( iom_use ('28d') ) THEN  ! depth of the 28 isotherm
             ztem2 = 28.
-            CALL dia_hth_dep( ztem2, hd28 )  
+            CALL dia_hth_dep( Kmm, ztem2, hd28 )  
             CALL iom_put( '28d', hd28 )    
          ENDIF
@@ -271 +260 @@
          IF( iom_use ('hc300') ) THEN  
             zzdep = 300.
-            CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc3 )
+            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc3 )
             CALL iom_put( 'hc300', rau0_rcp * htc3 )  ! vertically integrated heat content (J/m2)
          ENDIF
@@ -280 +269 @@
          IF( iom_use ('hc700') ) THEN  
             zzdep = 700.
-            CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc7 )
+            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc7 )
             CALL iom_put( 'hc700', rau0_rcp * htc7 )  ! vertically integrated heat content (J/m2)
   
@@ -290 +279 @@
          IF( iom_use ('hc2000') ) THEN  
             zzdep = 2000.
-            CALL  dia_hth_htc( zzdep, tsn(:,:,:,jp_tem), htc20 )
+            CALL  dia_hth_htc( Kmm, zzdep, ts(:,:,:,jp_tem,Kmm), htc20 )
             CALL iom_put( 'hc2000', rau0_rcp * htc20 )  ! vertically integrated heat content (J/m2)  
          ENDIF
@@ -301 +290 @@
    END SUBROUTINE dia_hth
 
-   SUBROUTINE dia_hth_dep( ptem, pdept )
-      !
+   SUBROUTINE dia_hth_dep( Kmm, ptem, pdept )
+      !
+      INTEGER , INTENT(in) :: Kmm      ! ocean time level index
       REAL(wp), INTENT(in) :: ptem
       REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pdept     
@@ -314 +304 @@
       ! --------------------------------------- !
       iktem(:,:) = 1
-      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zztmp = tsn(ji,jj,jk,jp_tem)
-               IF( zztmp >= ptem )   iktem(ji,jj) = jk
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         zztmp = ts(ji,jj,jk,jp_tem,Kmm)
+         IF( zztmp >= ptem )   iktem(ji,jj) = jk
+      END_3D
 
       ! ------------------------------- !
       !  Depth of ptem isotherm         !
       ! ------------------------------- !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean bottom
-            !
-            iid = iktem(ji,jj)
-            IF( iid /= 1 ) THEN 
-                zztmp =     gdept_n(ji,jj,iid  )   &                     ! linear interpolation
-                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
-                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
-                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
-               pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
-            ELSE 
-               pdept(ji,jj) = 0._wp
-            ENDIF
-         END DO
-      END DO
+      DO_2D_11_11
+         !
+         zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the ocean bottom
+         !
+         iid = iktem(ji,jj)
+         IF( iid /= 1 ) THEN 
+             zztmp =     gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation
+               &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   &
+               &  * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm)                       )   &
+               &  / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) )
+            pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
+         ELSE 
+            pdept(ji,jj) = 0._wp
+         ENDIF
+      END_2D
       !
    END SUBROUTINE dia_hth_dep
 
 
-   SUBROUTINE dia_hth_htc( pdep, ptn, phtc )
-      !
+   SUBROUTINE dia_hth_htc( Kmm, pdep, pt, phtc )
+      !
+      INTEGER , INTENT(in) :: Kmm      ! ocean time level index
       REAL(wp), INTENT(in) :: pdep     ! depth over the heat content
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: ptn   
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: pt   
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phtc  
       !
@@ -361 +346 @@
       
       IF( .NOT. ln_linssh ) THEN   ;   zthick(:,:) = 0._wp       ;   phtc(:,:) = 0._wp                                   
-      ELSE                         ;   zthick(:,:) = sshn(:,:)   ;   phtc(:,:) = ptn(:,:,1) * sshn(:,:) * tmask(:,:,1)   
+      ELSE                         ;   zthick(:,:) = ssh(:,:,Kmm)   ;   phtc(:,:) = pt(:,:,1) * ssh(:,:,Kmm) * tmask(:,:,1)   
       ENDIF
       !
       ilevel(:,:) = 1
-      DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( ( gdept_n(ji,jj,jk) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
-                   ilevel(ji,jj) = jk
-                   zthick(ji,jj) = zthick(ji,jj) + e3t_n(ji,jj,jk)
-                   phtc  (ji,jj) = phtc  (ji,jj) + e3t_n(ji,jj,jk) * ptn(ji,jj,jk)
-               ENDIF
-            ENDDO
-         ENDDO
-      ENDDO
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ik = ilevel(ji,jj)
-            zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep
-            phtc(ji,jj)   = phtc(ji,jj) + ptn(ji,jj,ik+1) * MIN( e3t_n(ji,jj,ik+1), zthick(ji,jj) ) &
-                                                          * tmask(ji,jj,ik+1)
-         END DO
-      ENDDO
+      DO_3D_11_11( 2, jpkm1 )
+         IF( ( gdept(ji,jj,jk,Kmm) < pdep ) .AND. ( tmask(ji,jj,jk) == 1 ) ) THEN
+             ilevel(ji,jj) = jk
+             zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
+             phtc  (ji,jj) = phtc  (ji,jj) + e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk)
+         ENDIF
+      END_3D
+      !
+      DO_2D_11_11
+         ik = ilevel(ji,jj)
+         zthick(ji,jj) = pdep - zthick(ji,jj)   !   remaining thickness to reach depht pdep
+         phtc(ji,jj)   = phtc(ji,jj) + pt(ji,jj,ik+1) * MIN( e3t(ji,jj,ik+1,Kmm), zthick(ji,jj) ) &
+                                                       * tmask(ji,jj,ik+1)
+      END_2D
       !
       !

NEMO/trunk/src/OCE/DIA/diaptr.F90

@@ -46 +46 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hstr_ove, hstr_btr, hstr_vtr   !: heat Salt TRansports(overturn, baro, merional)
 
-   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F)
-   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation
+   LOGICAL , PUBLIC ::   l_diaptr        !: tracers  trend flag (set from namelist in trdini)
    INTEGER, PARAMETER, PUBLIC ::   nptr = 5  ! (glo, atl, pac, ind, ipc)
 
@@ -60 +59 @@
    REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:) :: p_fval2d
 
+   LOGICAL ::   ll_init = .TRUE.        !: tracers  trend flag (set from namelist in trdini)
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -69 +69 @@
 CONTAINS
 
-   SUBROUTINE dia_ptr( pvtr )
+   SUBROUTINE dia_ptr( kt, Kmm, pvtr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dia_ptr  ***
       !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in)           ::   kt     ! ocean time-step index     
+      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
       !
@@ -92 +94 @@
       !
       IF( ln_timing )   CALL timing_start('dia_ptr')
-      !
+
+      IF( kt == nit000 .AND. ll_init )   CALL dia_ptr_init
+      !
+      IF( .NOT. l_diaptr )   RETURN
+
       IF( PRESENT( pvtr ) ) THEN
          IF( iom_use( 'zomsf' ) ) THEN    ! effective MSF
@@ -111 +117 @@
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpi
-                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
-                     zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
-                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
-                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                  ENDDO
-               ENDDO
-             ENDDO
+            DO_3D_10_11( 1, jpkm1 )
+               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+               zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
+               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
+               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
+            END_3D
          ENDIF
          IF( iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) ) THEN
@@ -186 +188 @@
          zts(:,:,:,:) = 0._wp
          IF( iom_use( 'zotem' ) .OR. iom_use( 'zosal' ) .OR. iom_use( 'zosrf' )  ) THEN    ! i-mean i-k-surface 
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk)
-                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
-                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
-                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
-                  END DO
-               END DO
-            END DO
+            DO_3D_11_11( 1, jpkm1 )
+               zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
+               zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
+               zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
+               zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc
+            END_3D
             !
             DO jn = 1, nptr
@@ -280 +278 @@
          IF( iom_use( 'sopstvtr' ) .OR. iom_use( 'sophtvtr' ) ) THEN
             zts(:,:,:,:) = 0._wp
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpi
-                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
-                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
-                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                  ENDDO
-               ENDDO
-             ENDDO
+            DO_3D_10_11( 1, jpkm1 )
+               zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
+               zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
+               zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
+            END_3D
              CALL dia_ptr_hst( jp_tem, 'vtr', zts(:,:,:,jp_tem) )
              CALL dia_ptr_hst( jp_sal, 'vtr', zts(:,:,:,jp_sal) )
@@ -326 +320 @@
       !! ** Purpose :   Initialization, namelist read
       !!----------------------------------------------------------------------
-      INTEGER ::  inum, jn, ios, ierr           ! local integers
-      !!
-      NAMELIST/namptr/ ln_diaptr, ln_subbas
+      INTEGER ::  inum, jn           ! local integers
+      !!
       REAL(wp), DIMENSION(jpi,jpj) :: zmsk
       !!----------------------------------------------------------------------
 
-
-      REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
-      READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist' )
-
-      REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
-      READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' )
-      IF(lwm) WRITE ( numond, namptr )
-
+      l_diaptr = .FALSE.
+      IF(   iom_use( 'zomsf'    ) .OR. iom_use( 'zotem'    ) .OR. iom_use( 'zosal'    ) .OR.  &
+         &  iom_use( 'zosrf'    ) .OR. iom_use( 'sopstove' ) .OR. iom_use( 'sophtove' ) .OR.  &
+         &  iom_use( 'sopstbtr' ) .OR. iom_use( 'sophtbtr' ) .OR. iom_use( 'sophtadv' ) .OR.  &
+         &  iom_use( 'sopstadv' ) .OR. iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) .OR.  & 
+         &  iom_use( 'sophteiv' ) .OR. iom_use( 'sopsteiv' ) .OR. iom_use( 'sopstvtr' ) .OR.  &
+         &  iom_use( 'sophtvtr' ) .OR. iom_use( 'uocetr_vsum_cumul' ) )  l_diaptr  = .TRUE.
+
+ 
       IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
@@ -347 +339 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namptr : set ptr parameters'
-         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr
+         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      l_diaptr  = ', l_diaptr
       ENDIF
 
-      IF( ln_diaptr ) THEN  
+      IF( l_diaptr ) THEN  
          !
          IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' )
@@ -389 +381 @@
          hstr_vtr(:,:,:) = 0._wp           !
          !
+         ll_init = .FALSE.
+         !
       ENDIF 
       ! 
@@ -394 +388 @@
 
 
-   SUBROUTINE dia_ptr_hst( ktra, cptr, pva ) 
+   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx ) 
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ptr_hst ***
@@ -403 +397 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx   ! 3D input array of advection/diffusion
       INTEGER                                        :: jn    !
 
@@ -410 +404 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_adv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_adv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_adv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_adv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -423 +417 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_ldf(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_ldf(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_ldf(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_ldf(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -436 +430 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_eiv(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_eiv(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_eiv(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_eiv(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -449 +443 @@
          IF( ktra == jp_tem )  THEN
              DO jn = 1, nptr
-                hstr_vtr(:,jp_tem,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_vtr(:,jp_tem,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
          IF( ktra == jp_sal )  THEN
              DO jn = 1, nptr
-                hstr_vtr(:,jp_sal,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                hstr_vtr(:,jp_sal,jn) = ptr_sj( pvflx(:,:,:), btmsk(:,:,jn) )
              ENDDO
          ENDIF
@@ -486 +480 @@
 
 
-   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sj_3d  ***
@@ -492 +486 @@
       !! ** Purpose :   i-k sum computation of a j-flux array
       !!
-      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
-      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
-      !!
-      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pva   ! mask flux array at V-point
+      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
+      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
+      !!
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)  ::   pvflx  ! mask flux array at V-point
       REAL(wp), INTENT(in), DIMENSION(jpi,jpj)      ::   pmsk   ! Optional 2D basin mask
       !
@@ -509 +503 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_00_00( 1, jpkm1 )
+         p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
+      END_3D
 #if defined key_mpp_mpi
       CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl)
@@ -523 +513 @@
 
 
-   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sj_2d  ***
       !!
-      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
-      !!
-      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
-      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
-      !!
-      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point
+      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array
+      !!
+      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
+      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
+      !!
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pvflx  ! mask flux array at V-point
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pmsk   ! Optional 2D basin mask
       !
@@ -546 +536 @@
       ijpj = jpj
       p_fval(:) = 0._wp
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! Vector opt.
-            p_fval(jj) = p_fval(jj) + pva(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
-         END DO
-      END DO
+      DO_2D_00_00
+         p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * pmsk(ji,jj) * tmask_i(ji,jj)
+      END_2D
 #if defined key_mpp_mpi
       CALL mpp_sum( 'diaptr', p_fval, ijpj, ncomm_znl )
@@ -577 +565 @@
       p_fval(:,:) = 0._wp
       DO jc = 1, jpnj ! looping over all processors in j axis
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
+         DO_2D_00_00
+            p_fval(ji,jj) = p_fval(ji,jj-1) + pva(ji,jj) * tmask_i(ji,jj)
+         END_2D
          CALL lbc_lnk( 'diaptr', p_fval, 'U', -1. )
       END DO
@@ -595 +581 @@
       !! ** Purpose :   i-sum computation of an array
       !!
-      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
-      !!
-      !! ** Action  : - p_fval: i-mean poleward flux of pva
+      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk).
+      !!
+      !! ** Action  : - p_fval: i-sum of masked field
       !!----------------------------------------------------------------------
       !!
@@ -618 +604 @@
       p_fval(:,:) = 0._wp
       !
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! Vector opt.
-               p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
-            END DO
-         END DO
-      END DO
+      DO_3D_00_00( 1, jpkm1 )
+         p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) * tmask_i(ji,jj)
+      END_3D
       !
 #if defined key_mpp_mpi

NEMO/trunk/src/OCE/DIA/diawri.F90

@@ -26 +26 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers 
+   USE isf_oce
+   USE isfcpl
+   USE abl            ! abl variables in case ln_abl = .true.
    USE dom_oce        ! ocean space and time domain
    USE phycst         ! physical constants
@@ -56 +59 @@
    USE lib_mpp         ! MPP library
    USE timing          ! preformance summary
-   USE diurnal_bulk    ! diurnal warm layer
-   USE cool_skin       ! Cool skin
+   USE diu_bulk        ! diurnal warm layer
+   USE diu_coolskin    ! Cool skin
 
    IMPLICIT NONE
@@ -65 +68 @@
    PUBLIC   dia_wri_state
    PUBLIC   dia_wri_alloc           ! Called by nemogcm module
-
+#if ! defined key_iomput   
+   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.)
+#endif
    INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
    INTEGER ::          nb_T              , ndim_bT   ! grid_T file
@@ -71 +76 @@
    INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
    INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
+   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file   
    INTEGER ::   ndex(1)                              ! ???
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
    INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
 
    !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -96 +103 @@
 
    
-   SUBROUTINE dia_wri( kt )
+   SUBROUTINE dia_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_wri  ***
@@ -106 +113 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
       !!
       INTEGER ::   ji, jj, jk       ! dummy loop indices
@@ -119 +127 @@
       ! Output the initial state and forcings
       IF( ninist == 1 ) THEN                       
-         CALL dia_wri_state( 'output.init' )
+         CALL dia_wri_state( Kmm, 'output.init' )
          ninist = 0
       ENDIF
@@ -128 +136 @@
       CALL iom_put("e3v_0", e3v_0(:,:,:) )
       !
-      CALL iom_put( "e3t" , e3t_n(:,:,:) )
-      CALL iom_put( "e3u" , e3u_n(:,:,:) )
-      CALL iom_put( "e3v" , e3v_n(:,:,:) )
-      CALL iom_put( "e3w" , e3w_n(:,:,:) )
+      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) )
+      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) )
+      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) )
+      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) )
       IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
+         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
 
       IF( ll_wd ) THEN
-         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
+         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
       ELSE
-         CALL iom_put( "ssh" , sshn )              ! sea surface height
+         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
       ENDIF
 
       IF( iom_use("wetdep") )   &                  ! wet depth
-         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
+         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )
       
-      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
-      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
+      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature
+      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature
       IF ( iom_use("sbt") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbkt(ji,jj)
-               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
-            END DO
-         END DO
+         DO_2D_11_11
+            ikbot = mbkt(ji,jj)
+            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
+         END_2D
          CALL iom_put( "sbt", z2d )                ! bottom temperature
       ENDIF
       
-      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
-      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
+      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity
+      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity
       IF ( iom_use("sbs") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbkt(ji,jj)
-               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
-            END DO
-         END DO
+         DO_2D_11_11
+            ikbot = mbkt(ji,jj)
+            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
+         END_2D
          CALL iom_put( "sbs", z2d )                ! bottom salinity
       ENDIF
@@ -171 +175 @@
          zztmp = rau0 * 0.25
          z2d(:,:) = 0._wp
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
-                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
-                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
-                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
-               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
-               !
-            END DO
-         END DO
+         DO_2D_00_00
+            zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   &
+               &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   &
+               &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   &
+               &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2
+            z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
+            !
+         END_2D
          CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
          CALL iom_put( "taubot", z2d )           
       ENDIF
          
-      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
-      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
+      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current
+      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current
       IF ( iom_use("sbu") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbku(ji,jj)
-               z2d(ji,jj) = un(ji,jj,ikbot)
-            END DO
-         END DO
+         DO_2D_11_11
+            ikbot = mbku(ji,jj)
+            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
+         END_2D
          CALL iom_put( "sbu", z2d )                ! bottom i-current
       ENDIF
       
-      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
-      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
+      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current
+      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current
       IF ( iom_use("sbv") ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ikbot = mbkv(ji,jj)
-               z2d(ji,jj) = vn(ji,jj,ikbot)
-            END DO
-         END DO
+         DO_2D_11_11
+            ikbot = mbkv(ji,jj)
+            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
+         END_2D
          CALL iom_put( "sbv", z2d )                ! bottom j-current
       ENDIF
 
-      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
-      !
-      CALL iom_put( "woce", wn )                   ! vertical velocity
+      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
+      !
+      CALL iom_put( "woce", ww )                   ! vertical velocity
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          z2d(:,:) = rau0 * e1e2t(:,:)
          DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
+            z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
          END DO
          CALL iom_put( "w_masstr" , z3d )  
@@ -222 +220 @@
       ENDIF
       !
-      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
+      IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
 
       CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
@@ -232 +230 @@
 
       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
-         DO jj = 2, jpjm1                                    ! sst gradient
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp  = tsn(ji,jj,1,jp_tem)
-               zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
-               zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
-               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
-                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
-            END DO
-         END DO
+         DO_2D_00_00
+            zztmp  = ts(ji,jj,1,jp_tem,Kmm)
+            zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
+            zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
+            z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
+               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
+         END_2D
          CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
          CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
@@ -250 +246 @@
       IF( iom_use("heatc") ) THEN
          z2d(:,:)  = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
+         END_3D
          CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
       ENDIF
@@ -262 +254 @@
       IF( iom_use("saltc") ) THEN
          z2d(:,:)  = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
+         END_3D
          CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
@@ -274 +262 @@
       IF ( iom_use("eken") ) THEN
          z3d(:,:,jpk) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
-                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
-                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
-                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+            z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   &
+               &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   &
+               &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   &
+               &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   )
+         END_3D
          CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
          CALL iom_put( "eken", z3d )                 ! kinetic energy
       ENDIF
       !
-      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
+      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence
       !
       IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
@@ -295 +279 @@
          z2d(:,:) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
+            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
             z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
@@ -304 +288 @@
       IF( iom_use("u_heattr") ) THEN
          z2d(:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
+         END_3D
          CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
          CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
@@ -317 +297 @@
       IF( iom_use("u_salttr") ) THEN
          z2d(:,:) = 0.e0 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
+         END_3D
          CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
          CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
@@ -332 +308 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
+            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
          END DO
          CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
@@ -339 +315 @@
       IF( iom_use("v_heattr") ) THEN
          z2d(:,:) = 0.e0 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
+         END_3D
          CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
          CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
@@ -352 +324 @@
       IF( iom_use("v_salttr") ) THEN
          z2d(:,:) = 0._wp 
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
+         END_3D
          CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
          CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
@@ -365 +333 @@
       IF( iom_use("tosmint") ) THEN
          z2d(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm)
+         END_3D
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
          CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
@@ -377 +341 @@
       IF( iom_use("somint") ) THEN
          z2d(:,:)=0._wp
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
+         END_3D
          CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
          CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
@@ -390 +350 @@
       CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
       !
-          
-      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
+      
+      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging
 
       IF( ln_timing )   CALL timing_stop('dia_wri')
@@ -411 +371 @@
          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
          !
-      dia_wri_alloc = MAXVAL(ierr)
+     dia_wri_alloc = MAXVAL(ierr)
       CALL mpp_sum( 'diawri', dia_wri_alloc )
       !
    END FUNCTION dia_wri_alloc
+ 
+   INTEGER FUNCTION dia_wri_alloc_abl()
+      !!----------------------------------------------------------------------
+     ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl)
+      CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
+      !
+   END FUNCTION dia_wri_alloc_abl
 
    
-   SUBROUTINE dia_wri( kt )
+   SUBROUTINE dia_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_wri  ***
@@ -430 +397 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
       !
       LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
@@ -437 +405 @@
       INTEGER  ::   ierr                                     ! error code return from allocation
       INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
+      INTEGER  ::   ipka                                     ! ABL
       INTEGER  ::   jn, ierror                               ! local integers
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
@@ -442 +411 @@
       REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
       !!----------------------------------------------------------------------
       !
       IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
-         CALL dia_wri_state( 'output.init' )
+         CALL dia_wri_state( Kmm, 'output.init' )
          ninist = 0
       ENDIF
@@ -475 +445 @@
       ijmi = 1      ;      ijma = jpj
       ipk = jpk
+      IF(ln_abl) ipka = jpkam1
 
       ! define time axis
@@ -577 +548 @@
             &          "m", ipk, gdepw_1d, nz_W, "down" )
 
+         IF( ln_abl ) THEN 
+         ! Define the ABL grid FILE ( nid_A )
+            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
+            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
+            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
+               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
+               &          nit000-1, zjulian, rdt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
+            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept
+               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" )
+            !                                                            ! Index of ocean points
+         ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 
+         zw3d_abl(:,:,:) = 1._wp 
+         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume
+            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface
+         DEALLOCATE(zw3d_abl)
+         ENDIF
 
          ! Declare all the output fields as NETCDF variables
@@ -586 +573 @@
             &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
          IF(  .NOT.ln_linssh  ) THEN
-            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
+            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
+            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
+            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
          ENDIF
@@ -607 +594 @@
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          IF(  ln_linssh  ) THEN
-            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
+            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm)
             &                                                                  , "KgC/m2/s",  &  ! sosst_cd
             &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
+            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm)
             &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
             &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -626 +613 @@
          CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-!
+         !
+         IF( ln_abl ) THEN
+            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
+            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
+            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl
+               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
+            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh
+               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                 
+#if defined key_si3
+            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i
+               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )
+#endif
+            CALL histend( nid_A, snc4chunks=snc4set )
+         ENDIF
+         !
          IF( ln_icebergs ) THEN
             CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
@@ -686 +697 @@
 
          !                                                                                      !!! nid_U : 3D
-         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
+         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm)
             &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
          IF( ln_wave .AND. ln_sdw) THEN
@@ -699 +710 @@
 
          !                                                                                      !!! nid_V : 3D
-         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
+         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm)
             &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
          IF( ln_wave .AND. ln_sdw) THEN
@@ -712 +723 @@
 
          !                                                                                      !!! nid_W : 3D
-         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
+         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww
             &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
          CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
@@ -750 +761 @@
 
       IF( .NOT.ln_linssh ) THEN
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
+         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content
       ELSE
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
+         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
+         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
+         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
+         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
       ENDIF
       IF( .NOT.ln_linssh ) THEN
-         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
-         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
+         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness
+         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
          CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
-      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
+      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
       CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
       CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
@@ -773 +784 @@
                                                                                   ! in linear free surface case)
       IF( ln_linssh ) THEN
-         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
+         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
          CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
-         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
+         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
          CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
       ENDIF
@@ -784 +795 @@
       CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
       CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
-!
+      !
+      IF( ln_abl ) THEN 
+         ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
+         IF( ln_mskland )   THEN 
+            DO jk=1,jpka
+               zw3d_abl(:,:,jk) = tmask(:,:,1)
+            END DO       
+         ELSE
+            zw3d_abl(:,:,:) = 1._wp     
+         ENDIF       
+         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh 
+         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
+         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
+         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
+         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl       
+         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
+         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
+         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl 
+#if defined key_si3
+         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
+#endif
+         DEALLOCATE(zw3d_abl)
+      ENDIF
+      !
       IF( ln_icebergs ) THEN
          !
@@ -811 +845 @@
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+         zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
@@ -824 +858 @@
 #endif
 
-      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
+      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
       CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
-      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
+      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
       IF( ln_zad_Aimp ) THEN
-         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current
+         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current
       ELSE
-         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current
+         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
       ENDIF
       CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
@@ -854 +888 @@
          CALL histclo( nid_V )
          CALL histclo( nid_W )
+         IF(ln_abl) CALL histclo( nid_A )
       ENDIF
       !
@@ -861 +896 @@
 #endif
 
-   SUBROUTINE dia_wri_state( cdfile_name )
+   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE dia_wri_state  ***
@@ -874 +909 @@
       !!      File 'output.abort.nc' is created in case of abnormal job end
       !!----------------------------------------------------------------------
+      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
       CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
       !!
-      INTEGER :: inum
+      INTEGER :: inum, jk
       !!----------------------------------------------------------------------
       ! 
@@ -890 +926 @@
 #endif
 
-      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
-      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
-      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
-      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
-      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
+      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
+      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
+      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
+      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
+      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
       IF( ln_zad_Aimp ) THEN
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
       ELSE
-         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity
-      ENDIF
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
+      ENDIF
+      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
+      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height
+
+      IF ( ln_isf ) THEN
+         IF (ln_isfcav_mlt) THEN
+            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 )
+         END IF
+         IF (ln_isfpar_mlt) THEN
+            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,8)  )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 )
+         END IF
+      END IF
+
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
@@ -915 +974 @@
       CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
       IF(  .NOT.ln_linssh  ) THEN             
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth 
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness  
+         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth 
+         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness  
       END IF
       IF( ln_wave .AND. ln_sdw ) THEN
@@ -923 +982 @@
          CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
       ENDIF
+      IF ( ln_abl ) THEN
+         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind
+         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind
+         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature
+         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity
+      ENDIF
  
 #if defined key_si3

@@ -3 +3 @@
 ^/utils/build/mk@HEAD         mk
 ^/utils/tools@HEAD            tools
-^/vendors/AGRIF/dev@HEAD      ext/AGRIF
+^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF
 ^/vendors/FCM@HEAD            ext/FCM
 ^/vendors/IOIPSL@HEAD         ext/IOIPSL