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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12928 for NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90 – NEMO

Ignore:
Timestamp:
2020-05-14T21:46:00+02:00 (4 years ago)
Author:
smueller
Message:

Synchronizing with /NEMO/trunk@12925 (ticket #2170)

Location:
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser

    • Property svn:externals
      •  

        old new  
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@HEAD         sette 
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90

    r12383 r12928  
    2626   !!---------------------------------------------------------------------- 
    2727   USE oce            ! ocean dynamics and tracers  
     28   USE isf_oce 
     29   USE isfcpl 
     30   USE abl            ! abl variables in case ln_abl = .true. 
    2831   USE dom_oce        ! ocean space and time domain 
    2932   USE phycst         ! physical constants 
     
    4750   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    4851   USE in_out_manager ! I/O manager 
    49    USE diatmb         ! Top,middle,bottom output 
    5052   USE dia25h         ! 25h Mean output 
    5153   USE iom            !  
     
    5860   USE lib_mpp         ! MPP library 
    5961   USE timing          ! preformance summary 
    60    USE diurnal_bulk    ! diurnal warm layer 
    61    USE cool_skin       ! Cool skin 
     62   USE diu_bulk        ! diurnal warm layer 
     63   USE diu_coolskin    ! Cool skin 
    6264 
    6365   IMPLICIT NONE 
     
    6769   PUBLIC   dia_wri_state 
    6870   PUBLIC   dia_wri_alloc           ! Called by nemogcm module 
    69  
     71#if ! defined key_iomput    
     72   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.) 
     73#endif 
    7074   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
    7175   INTEGER ::          nb_T              , ndim_bT   ! grid_T file 
     
    7377   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
    7478   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
     79   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7580   INTEGER ::   ndex(1)                              ! ??? 
    7681   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     82   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 
    7783   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
    7884   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    7985 
    8086   !! * Substitutions 
    81 #  include "vectopt_loop_substitute.h90" 
     87#  include "do_loop_substitute.h90" 
    8288   !!---------------------------------------------------------------------- 
    8389   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    98104 
    99105    
    100    SUBROUTINE dia_wri( kt ) 
     106   SUBROUTINE dia_wri( kt, Kmm ) 
    101107      !!--------------------------------------------------------------------- 
    102108      !!                  ***  ROUTINE dia_wri  *** 
     
    108114      !!---------------------------------------------------------------------- 
    109115      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     116      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index 
    110117      !! 
    111118      INTEGER ::   ji, jj, jk       ! dummy loop indices 
     
    121128      ! Output the initial state and forcings 
    122129      IF( ninist == 1 ) THEN                        
    123          CALL dia_wri_state( 'output.init' ) 
     130         CALL dia_wri_state( Kmm, 'output.init' ) 
    124131         ninist = 0 
    125132      ENDIF 
     
    130137      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
    131138      ! 
    132       CALL iom_put( "e3t" , e3t_n(:,:,:) ) 
    133       CALL iom_put( "e3u" , e3u_n(:,:,:) ) 
    134       CALL iom_put( "e3v" , e3v_n(:,:,:) ) 
    135       CALL iom_put( "e3w" , e3w_n(:,:,:) ) 
     139      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 
     140      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 
     141      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 
     142      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 
    136143      IF( iom_use("e3tdef") )   & 
    137          CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
     144         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    138145 
    139146      IF( ll_wd ) THEN 
    140          CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
     147         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
    141148      ELSE 
    142          CALL iom_put( "ssh" , sshn )              ! sea surface height 
     149         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    143150      ENDIF 
    144151 
    145152      IF( iom_use("wetdep") )   &                  ! wet depth 
    146          CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 
     153         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 
    147154       
    148       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
    149       CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
     155      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature 
     156      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    150157      IF ( iom_use("sbt") ) THEN 
    151          DO jj = 1, jpj 
    152             DO ji = 1, jpi 
    153                ikbot = mbkt(ji,jj) 
    154                z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem) 
    155             END DO 
    156          END DO 
     158         DO_2D_11_11 
     159            ikbot = mbkt(ji,jj) 
     160            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
     161         END_2D 
    157162         CALL iom_put( "sbt", z2d )                ! bottom temperature 
    158163      ENDIF 
    159164       
    160       CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
    161       CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
     165      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity 
     166      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    162167      IF ( iom_use("sbs") ) THEN 
    163          DO jj = 1, jpj 
    164             DO ji = 1, jpi 
    165                ikbot = mbkt(ji,jj) 
    166                z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal) 
    167             END DO 
    168          END DO 
     168         DO_2D_11_11 
     169            ikbot = mbkt(ji,jj) 
     170            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
     171         END_2D 
    169172         CALL iom_put( "sbs", z2d )                ! bottom salinity 
    170173      ENDIF 
    171174 
    172175      IF ( iom_use("taubot") ) THEN                ! bottom stress 
    173          zztmp = rau0 * 0.25 
     176         zztmp = rho0 * 0.25 
    174177         z2d(:,:) = 0._wp 
    175          DO jj = 2, jpjm1 
    176             DO ji = fs_2, fs_jpim1   ! vector opt. 
    177                zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   & 
    178                   &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   & 
    179                   &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   & 
    180                   &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2 
    181                z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
    182                ! 
    183             END DO 
    184          END DO 
     178         DO_2D_00_00 
     179            zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
     180               &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
     181               &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   & 
     182               &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2 
     183            z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
     184            ! 
     185         END_2D 
    185186         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    186187         CALL iom_put( "taubot", z2d )            
    187188      ENDIF 
    188189          
    189       CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current 
    190       CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current 
     190      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current 
     191      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    191192      IF ( iom_use("sbu") ) THEN 
    192          DO jj = 1, jpj 
    193             DO ji = 1, jpi 
    194                ikbot = mbku(ji,jj) 
    195                z2d(ji,jj) = un(ji,jj,ikbot) 
    196             END DO 
    197          END DO 
     193         DO_2D_11_11 
     194            ikbot = mbku(ji,jj) 
     195            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     196         END_2D 
    198197         CALL iom_put( "sbu", z2d )                ! bottom i-current 
    199198      ENDIF 
    200199       
    201       CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current 
    202       CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current 
     200      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current 
     201      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    203202      IF ( iom_use("sbv") ) THEN 
    204          DO jj = 1, jpj 
    205             DO ji = 1, jpi 
    206                ikbot = mbkv(ji,jj) 
    207                z2d(ji,jj) = vn(ji,jj,ikbot) 
    208             END DO 
    209          END DO 
     203         DO_2D_11_11 
     204            ikbot = mbkv(ji,jj) 
     205            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
     206         END_2D 
    210207         CALL iom_put( "sbv", z2d )                ! bottom j-current 
    211208      ENDIF 
    212209 
    213       IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
    214       ! 
    215       CALL iom_put( "woce", wn )                   ! vertical velocity 
     210      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
     211      ! 
     212      CALL iom_put( "woce", ww )                   ! vertical velocity 
    216213      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
    217214         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    218          z2d(:,:) = rau0 * e1e2t(:,:) 
     215         z2d(:,:) = rho0 * e1e2t(:,:) 
    219216         DO jk = 1, jpk 
    220             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
     217            z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:) 
    221218         END DO 
    222219         CALL iom_put( "w_masstr" , z3d )   
     
    224221      ENDIF 
    225222      ! 
    226       IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
     223      IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
    227224 
    228225      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     
    234231 
    235232      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    236          DO jj = 2, jpjm1                                    ! sst gradient 
    237             DO ji = fs_2, fs_jpim1   ! vector opt. 
    238                zztmp  = tsn(ji,jj,1,jp_tem) 
    239                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) 
    240                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) 
    241                z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    242                   &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    243             END DO 
    244          END DO 
     233         DO_2D_00_00 
     234            zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
     235            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) 
     236            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) 
     237            z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
     238               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
     239         END_2D 
    245240         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    246241         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
     
    252247      IF( iom_use("heatc") ) THEN 
    253248         z2d(:,:)  = 0._wp  
    254          DO jk = 1, jpkm1 
    255             DO jj = 1, jpj 
    256                DO ji = 1, jpi 
    257                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
    258                END DO 
    259             END DO 
    260          END DO 
    261          CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
     249         DO_3D_11_11( 1, jpkm1 ) 
     250            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
     251         END_3D 
     252         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    262253      ENDIF 
    263254 
    264255      IF( iom_use("saltc") ) THEN 
    265256         z2d(:,:)  = 0._wp  
    266          DO jk = 1, jpkm1 
    267             DO jj = 1, jpj 
    268                DO ji = 1, jpi 
    269                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
    270                END DO 
    271             END DO 
    272          END DO 
    273          CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
     257         DO_3D_11_11( 1, jpkm1 ) 
     258            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
     259         END_3D 
     260         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    274261      ENDIF 
    275262      ! 
    276263      IF ( iom_use("eken") ) THEN 
    277264         z3d(:,:,jpk) = 0._wp  
    278          DO jk = 1, jpkm1 
    279             DO jj = 2, jpjm1 
    280                DO ji = fs_2, fs_jpim1   ! vector opt. 
    281                   zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    282                   z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
    283                      &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
    284                      &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
    285                      &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
    286                END DO 
    287             END DO 
    288          END DO 
     265         DO_3D_00_00( 1, jpkm1 ) 
     266            zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     267            z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
     268               &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
     269               &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
     270               &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
     271         END_3D 
    289272         CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    290273         CALL iom_put( "eken", z3d )                 ! kinetic energy 
    291274      ENDIF 
    292275      ! 
    293       CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence 
     276      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
    294277      ! 
    295278      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
     
    297280         z2d(:,:) = 0.e0 
    298281         DO jk = 1, jpkm1 
    299             z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 
     282            z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk) 
    300283            z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    301284         END DO 
     
    306289      IF( iom_use("u_heattr") ) THEN 
    307290         z2d(:,:) = 0._wp  
    308          DO jk = 1, jpkm1 
    309             DO jj = 2, jpjm1 
    310                DO ji = fs_2, fs_jpim1   ! vector opt. 
    311                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
    312                END DO 
    313             END DO 
    314          END DO 
     291         DO_3D_00_00( 1, jpkm1 ) 
     292            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) ) 
     293         END_3D 
    315294         CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    316295         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
     
    319298      IF( iom_use("u_salttr") ) THEN 
    320299         z2d(:,:) = 0.e0  
    321          DO jk = 1, jpkm1 
    322             DO jj = 2, jpjm1 
    323                DO ji = fs_2, fs_jpim1   ! vector opt. 
    324                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
    325                END DO 
    326             END DO 
    327          END DO 
     300         DO_3D_00_00( 1, jpkm1 ) 
     301            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) ) 
     302         END_3D 
    328303         CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    329304         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
     
    334309         z3d(:,:,jpk) = 0.e0 
    335310         DO jk = 1, jpkm1 
    336             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 
     311            z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    337312         END DO 
    338313         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction 
     
    341316      IF( iom_use("v_heattr") ) THEN 
    342317         z2d(:,:) = 0.e0  
    343          DO jk = 1, jpkm1 
    344             DO jj = 2, jpjm1 
    345                DO ji = fs_2, fs_jpim1   ! vector opt. 
    346                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
    347                END DO 
    348             END DO 
    349          END DO 
     318         DO_3D_00_00( 1, jpkm1 ) 
     319            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) ) 
     320         END_3D 
    350321         CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    351322         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
     
    354325      IF( iom_use("v_salttr") ) THEN 
    355326         z2d(:,:) = 0._wp  
    356          DO jk = 1, jpkm1 
    357             DO jj = 2, jpjm1 
    358                DO ji = fs_2, fs_jpim1   ! vector opt. 
    359                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
    360                END DO 
    361             END DO 
    362          END DO 
     327         DO_3D_00_00( 1, jpkm1 ) 
     328            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) ) 
     329         END_3D 
    363330         CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    364331         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
     
    367334      IF( iom_use("tosmint") ) THEN 
    368335         z2d(:,:) = 0._wp 
    369          DO jk = 1, jpkm1 
    370             DO jj = 2, jpjm1 
    371                DO ji = fs_2, fs_jpim1   ! vector opt. 
    372                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
    373                END DO 
    374             END DO 
    375          END DO 
     336         DO_3D_00_00( 1, jpkm1 ) 
     337            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
     338         END_3D 
    376339         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    377          CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature 
     340         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
    378341      ENDIF 
    379342      IF( iom_use("somint") ) THEN 
    380343         z2d(:,:)=0._wp 
    381          DO jk = 1, jpkm1 
    382             DO jj = 2, jpjm1 
    383                DO ji = fs_2, fs_jpim1   ! vector opt. 
    384                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
    385                END DO 
    386             END DO 
    387          END DO 
     344         DO_3D_00_00( 1, jpkm1 ) 
     345            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
     346         END_3D 
    388347         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    389          CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity 
     348         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
    390349      ENDIF 
    391350 
    392351      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2) 
    393352      ! 
    394  
    395       IF (ln_diatmb)   CALL dia_tmb                   ! tmb values  
    396            
    397       IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging 
     353       
     354      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
    398355 
    399356      IF( ln_timing )   CALL timing_stop('dia_wri') 
     
    410367      INTEGER, DIMENSION(2) :: ierr 
    411368      !!---------------------------------------------------------------------- 
    412       ierr = 0 
    413       ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
    414          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    415          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
     369      IF( nn_write == -1 ) THEN 
     370         dia_wri_alloc = 0 
     371      ELSE     
     372         ierr = 0 
     373         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
     374            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
     375            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
    416376         ! 
    417       dia_wri_alloc = MAXVAL(ierr) 
    418       CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     377         dia_wri_alloc = MAXVAL(ierr) 
     378         CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     379         ! 
     380      ENDIF 
    419381      ! 
    420382   END FUNCTION dia_wri_alloc 
     383  
     384   INTEGER FUNCTION dia_wri_alloc_abl() 
     385      !!---------------------------------------------------------------------- 
     386     ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl) 
     387      CALL mpp_sum( 'diawri', dia_wri_alloc_abl ) 
     388      ! 
     389   END FUNCTION dia_wri_alloc_abl 
    421390 
    422391    
    423    SUBROUTINE dia_wri( kt ) 
     392   SUBROUTINE dia_wri( kt, Kmm ) 
    424393      !!--------------------------------------------------------------------- 
    425394      !!                  ***  ROUTINE dia_wri  *** 
     
    434403      !!---------------------------------------------------------------------- 
    435404      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     405      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index 
    436406      ! 
    437407      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout 
     
    441411      INTEGER  ::   ierr                                     ! error code return from allocation 
    442412      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers 
     413      INTEGER  ::   ipka                                     ! ABL 
    443414      INTEGER  ::   jn, ierror                               ! local integers 
    444415      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars 
     
    446417      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    447418      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     419      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace 
    448420      !!---------------------------------------------------------------------- 
    449421      ! 
    450422      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==! 
    451          CALL dia_wri_state( 'output.init' ) 
     423         CALL dia_wri_state( Kmm, 'output.init' ) 
    452424         ninist = 0 
    453425      ENDIF 
     
    466438      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes) 
    467439#if defined key_diainstant 
    468       zsto = nn_write * rdt 
     440      zsto = nn_write * rn_Dt 
    469441      clop = "inst("//TRIM(clop)//")" 
    470442#else 
    471       zsto=rdt 
     443      zsto=rn_Dt 
    472444      clop = "ave("//TRIM(clop)//")" 
    473445#endif 
    474       zout = nn_write * rdt 
    475       zmax = ( nitend - nit000 + 1 ) * rdt 
     446      zout = nn_write * rn_Dt 
     447      zmax = ( nitend - nit000 + 1 ) * rn_Dt 
    476448 
    477449      ! Define indices of the horizontal output zoom and vertical limit storage 
     
    479451      ijmi = 1      ;      ijma = jpj 
    480452      ipk = jpk 
     453      IF(ln_abl) ipka = jpkam1 
    481454 
    482455      ! define time axis 
     
    493466 
    494467         ! Compute julian date from starting date of the run 
    495          CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) 
     468         CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian ) 
    496469         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment 
    497470         IF(lwp)WRITE(numout,*) 
     
    515488         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
    516489            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    517             &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 
     490            &          nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 
    518491         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept 
    519492            &           "m", ipk, gdept_1d, nz_T, "down" ) 
     
    551524         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu 
    552525            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    553             &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 
     526            &          nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 
    554527         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept 
    555528            &           "m", ipk, gdept_1d, nz_U, "down" ) 
     
    564537         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv 
    565538            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    566             &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 
     539            &          nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 
    567540         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept 
    568541            &          "m", ipk, gdept_1d, nz_V, "down" ) 
     
    577550         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
    578551            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    579             &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 
     552            &          nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 
    580553         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw 
    581554            &          "m", ipk, gdepw_1d, nz_W, "down" ) 
    582555 
     556         IF( ln_abl ) THEN  
     557         ! Define the ABL grid FILE ( nid_A ) 
     558            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' ) 
     559            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
     560            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     561               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     562               &          nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set ) 
     563            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept 
     564               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" ) 
     565            !                                                            ! Index of ocean points 
     566         ALLOCATE( zw3d_abl(jpi,jpj,ipka) )  
     567         zw3d_abl(:,:,:) = 1._wp  
     568         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume 
     569            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface 
     570         DEALLOCATE(zw3d_abl) 
     571         ENDIF 
    583572 
    584573         ! Declare all the output fields as NETCDF variables 
     
    590579            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    591580         IF(  .NOT.ln_linssh  ) THEN 
    592             CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n 
     581            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm) 
    593582            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    594             CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n 
     583            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm) 
    595584            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    596             CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n 
     585            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm) 
    597586            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    598587         ENDIF 
     
    611600            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    612601         IF(  ln_linssh  ) THEN 
    613             CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem) 
     602            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm) 
    614603            &                                                                  , "KgC/m2/s",  &  ! sosst_cd 
    615604            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    616             CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal) 
     605            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm) 
    617606            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd 
    618607            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    630619         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm 
    631620            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    632 ! 
     621         ! 
     622         IF( ln_abl ) THEN 
     623            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl 
     624               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     625            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl 
     626               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     627            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl 
     628               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     629            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl 
     630               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     631            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl 
     632               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     633            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl 
     634               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     635            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl 
     636               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     637            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh 
     638               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                  
     639#if defined key_si3 
     640            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i 
     641               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout ) 
     642#endif 
     643            CALL histend( nid_A, snc4chunks=snc4set ) 
     644         ENDIF 
     645         ! 
    633646         IF( ln_icebergs ) THEN 
    634647            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , & 
     
    690703 
    691704         !                                                                                      !!! nid_U : 3D 
    692          CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un 
     705         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm) 
    693706            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
    694707         IF( ln_wave .AND. ln_sdw) THEN 
     
    703716 
    704717         !                                                                                      !!! nid_V : 3D 
    705          CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn 
     718         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm) 
    706719            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
    707720         IF( ln_wave .AND. ln_sdw) THEN 
     
    716729 
    717730         !                                                                                      !!! nid_W : 3D 
    718          CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn 
     731         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww 
    719732            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    720733         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt 
     
    754767 
    755768      IF( .NOT.ln_linssh ) THEN 
    756          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content 
    757          CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content 
    758          CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content 
    759          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content 
     769         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content 
     770         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content 
     771         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content 
     772         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content 
    760773      ELSE 
    761          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature 
    762          CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity 
    763          CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature 
    764          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity 
     774         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature 
     775         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity 
     776         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature 
     777         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity 
    765778      ENDIF 
    766779      IF( .NOT.ln_linssh ) THEN 
    767          zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    768          CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness 
    769          CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth 
     780         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
     781         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness 
     782         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth 
    770783         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation 
    771784      ENDIF 
    772       CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height 
     785      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height 
    773786      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux 
    774787      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs 
     
    777790                                                                                  ! in linear free surface case) 
    778791      IF( ln_linssh ) THEN 
    779          zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem) 
     792         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm) 
    780793         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst 
    781          zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal) 
     794         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm) 
    782795         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss 
    783796      ENDIF 
     
    788801      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction    
    789802      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed    
    790 ! 
     803      ! 
     804      IF( ln_abl ) THEN  
     805         ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 
     806         IF( ln_mskland )   THEN  
     807            DO jk=1,jpka 
     808               zw3d_abl(:,:,jk) = tmask(:,:,1) 
     809            END DO        
     810         ELSE 
     811            zw3d_abl(:,:,:) = 1._wp      
     812         ENDIF        
     813         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh  
     814         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl 
     815         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl 
     816         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl 
     817         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl        
     818         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl 
     819         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl 
     820         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl  
     821#if defined key_si3 
     822         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i 
     823#endif 
     824         DEALLOCATE(zw3d_abl) 
     825      ENDIF 
     826      ! 
    791827      IF( ln_icebergs ) THEN 
    792828         ! 
     
    815851         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    816852         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    817          zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     853         zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 
    818854         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    819855      ENDIF 
     
    828864#endif 
    829865 
    830       CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current 
     866      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current 
    831867      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress 
    832868 
    833       CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current 
     869      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current 
    834870      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    835871 
    836872      IF( ln_zad_Aimp ) THEN 
    837          CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current 
     873         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current 
    838874      ELSE 
    839          CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current 
     875         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current 
    840876      ENDIF 
    841877      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
     
    858894         CALL histclo( nid_V ) 
    859895         CALL histclo( nid_W ) 
     896         IF(ln_abl) CALL histclo( nid_A ) 
    860897      ENDIF 
    861898      ! 
     
    865902#endif 
    866903 
    867    SUBROUTINE dia_wri_state( cdfile_name ) 
     904   SUBROUTINE dia_wri_state( Kmm, cdfile_name ) 
    868905      !!--------------------------------------------------------------------- 
    869906      !!                 ***  ROUTINE dia_wri_state  *** 
     
    878915      !!      File 'output.abort.nc' is created in case of abnormal job end 
    879916      !!---------------------------------------------------------------------- 
     917      INTEGER           , INTENT( in ) ::   Kmm              ! time level index 
    880918      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created 
    881919      !! 
    882       INTEGER :: inum 
     920      INTEGER :: inum, jk 
    883921      !!---------------------------------------------------------------------- 
    884922      !  
     
    887925      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
    888926      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
    889  
    890 #if defined key_si3 
    891      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
    892 #else 
    893      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
    894 #endif 
    895  
    896       CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature 
    897       CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity 
    898       CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height 
    899       CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity 
    900       CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity 
     927      ! 
     928      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
     929      ! 
     930      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature 
     931      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity 
     932      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height 
     933      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity 
     934      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity 
    901935      IF( ln_zad_Aimp ) THEN 
    902          CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity 
     936         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity 
    903937      ELSE 
    904          CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity 
    905       ENDIF 
     938         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity 
     939      ENDIF 
     940      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity 
     941      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height 
     942      ! 
     943      IF ( ln_isf ) THEN 
     944         IF (ln_isfcav_mlt) THEN 
     945            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity 
     946            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity 
     947            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity 
     948            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity 
     949            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity 
     950            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 
     951         END IF 
     952         IF (ln_isfpar_mlt) THEN 
     953            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity 
     954            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity 
     955            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity 
     956            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity 
     957            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity 
     958            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity 
     959            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 
     960         END IF 
     961      END IF 
     962      ! 
    906963      IF( ALLOCATED(ahtu) ) THEN 
    907964         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
     
    919976      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress 
    920977      IF(  .NOT.ln_linssh  ) THEN              
    921          CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth  
    922          CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness   
     978         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth  
     979         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness   
    923980      END IF 
    924981      IF( ln_wave .AND. ln_sdw ) THEN 
     
    942999      ENDIF 
    9431000  
     1001      IF ( ln_abl ) THEN 
     1002         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind 
     1003         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind 
     1004         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature 
     1005         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity 
     1006      ENDIF 
     1007      ! 
     1008      CALL iom_close( inum ) 
     1009      !  
    9441010#if defined key_si3 
    9451011      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid 
     1012         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 
    9461013         CALL ice_wri_state( inum ) 
     1014         CALL iom_close( inum ) 
    9471015      ENDIF 
    9481016#endif 
    949       ! 
    950       CALL iom_close( inum ) 
    951       !  
     1017 
    9521018   END SUBROUTINE dia_wri_state 
    9531019 
Note: See TracChangeset for help on using the changeset viewer.