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 12740 for NEMO/trunk/tests/CANAL/MY_SRC/diawri.F90 – NEMO

Ignore:
Timestamp:
2020-04-12T11:03:06+02:00 (4 years ago)
Author:
smasson
Message:

trunk: update/debug of tests and C1D, see #2442

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/tests/CANAL/MY_SRC/diawri.F90

    r12489 r12740  
    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 
     
    6568   PUBLIC   dia_wri_state 
    6669   PUBLIC   dia_wri_alloc           ! Called by nemogcm module 
    67  
     70#if ! defined key_iomput    
     71   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.) 
     72#endif 
    6873   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file 
    6974   INTEGER ::          nb_T              , ndim_bT   ! grid_T file 
     
    7176   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file 
    7277   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file 
     78   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file    
    7379   INTEGER ::   ndex(1)                              ! ??? 
    7480   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 
     81   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL 
    7582   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 
    7683   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT 
    7784 
     85   !! * Substitutions 
     86#  include "do_loop_substitute.h90" 
    7887   !!---------------------------------------------------------------------- 
    7988   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    147156      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    148157      IF ( iom_use("sbt") ) THEN 
    149          DO jj = 1, jpj 
    150             DO ji = 1, jpi 
    151                ikbot = mbkt(ji,jj) 
    152                z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
    153             END DO 
    154          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 
    155162         CALL iom_put( "sbt", z2d )                ! bottom temperature 
    156163      ENDIF 
     
    159166      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    160167      IF ( iom_use("sbs") ) THEN 
    161          DO jj = 1, jpj 
    162             DO ji = 1, jpi 
    163                ikbot = mbkt(ji,jj) 
    164                z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
    165             END DO 
    166          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 
    167172         CALL iom_put( "sbs", z2d )                ! bottom salinity 
    168173      ENDIF 
     
    171176         zztmp = rho0 * 0.25 
    172177         z2d(:,:) = 0._wp 
    173          DO jj = 2, jpjm1 
    174             DO ji = fs_2, fs_jpim1   ! vector opt. 
    175                zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
    176                   &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
    177                   &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   & 
    178                   &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2 
    179                z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
    180                ! 
    181             END DO 
    182          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 
    183186         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    184187         CALL iom_put( "taubot", z2d )            
     
    188191      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    189192      IF ( iom_use("sbu") ) THEN 
    190          DO jj = 1, jpj 
    191             DO ji = 1, jpi 
    192                ikbot = mbku(ji,jj) 
    193                z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
    194             END DO 
    195          END DO 
     193         DO_2D_11_11 
     194            ikbot = mbku(ji,jj) 
     195            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
     196         END_2D 
    196197         CALL iom_put( "sbu", z2d )                ! bottom i-current 
    197198      ENDIF 
     
    200201      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    201202      IF ( iom_use("sbv") ) THEN 
    202          DO jj = 1, jpj 
    203             DO ji = 1, jpi 
    204                ikbot = mbkv(ji,jj) 
    205                z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
    206             END DO 
    207          END DO 
     203         DO_2D_11_11 
     204            ikbot = mbkv(ji,jj) 
     205            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
     206         END_2D 
    208207         CALL iom_put( "sbv", z2d )                ! bottom j-current 
    209208      ENDIF 
    210209 
     210      IF( ln_zad_Aimp ) ww = ww + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output 
     211      ! 
    211212      CALL iom_put( "woce", ww )                   ! vertical velocity 
    212213      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
     
    219220         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
    220221      ENDIF 
     222      ! 
     223      IF( ln_zad_Aimp ) ww = ww - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output 
    221224 
    222225      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     
    227230      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    228231 
    229       IF ( iom_use("salgrad") .OR. iom_use("salgrad2") ) THEN 
    230          z3d(:,:,jpk) = 0. 
    231          DO jk = 1, jpkm1 
    232             DO jj = 2, jpjm1                                    ! sal gradient 
    233                DO ji = fs_2, fs_jpim1   ! vector opt. 
    234                   zztmp  = ts(ji,jj,jk,jp_sal,Kmm) 
    235                   zztmpx = ( ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm) ) * r1_e1u(ji-1,jj) 
    236                   zztmpy = ( ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm) ) * r1_e2v(ji,jj-1) 
    237                   z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    238                      &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk) 
    239                END DO 
    240             END DO 
    241          END DO 
    242          CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    243          CALL iom_put( "salgrad2",  z3d )          ! square of module of sal gradient 
    244          z3d(:,:,:) = SQRT( z3d(:,:,:) ) 
    245          CALL iom_put( "salgrad" ,  z3d )          ! module of sal gradient 
    246       ENDIF 
    247           
    248232      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    249          DO jj = 2, jpjm1                                    ! sst gradient 
    250             DO ji = fs_2, fs_jpim1   ! vector opt. 
    251                zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
    252                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) 
    253                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) 
    254                z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    255                   &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
    256             END DO 
    257          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 
    258240         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    259241         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
     
    265247      IF( iom_use("heatc") ) THEN 
    266248         z2d(:,:)  = 0._wp  
    267          DO jk = 1, jpkm1 
    268             DO jj = 1, jpj 
    269                DO ji = 1, jpi 
    270                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
    271                END DO 
    272             END DO 
    273          END DO 
     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 
    274252         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    275253      ENDIF 
     
    277255      IF( iom_use("saltc") ) THEN 
    278256         z2d(:,:)  = 0._wp  
    279          DO jk = 1, jpkm1 
    280             DO jj = 1, jpj 
    281                DO ji = 1, jpi 
    282                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    283                END DO 
    284             END DO 
    285          END DO 
     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 
    286260         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    287261      ENDIF 
     
    289263      IF( iom_use("salt2c") ) THEN 
    290264         z2d(:,:)  = 0._wp  
    291          DO jk = 1, jpkm1 
    292             DO jj = 1, jpj 
    293                DO ji = 1, jpi 
    294                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    295                END DO 
    296             END DO 
    297          END DO 
     265         DO_3D_11_11( 1, jpkm1 ) 
     266            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
     267         END_3D 
    298268         CALL iom_put( "salt2c", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    299269      ENDIF 
     
    301271      IF ( iom_use("eken") ) THEN 
    302272         z3d(:,:,jpk) = 0._wp  
    303          DO jk = 1, jpkm1 
    304             DO jj = 2, jpjm1 
    305                DO ji = fs_2, fs_jpim1   ! vector opt. 
    306                   zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    307                   z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
    308                      &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
    309                      &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
    310                      &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
    311                END DO 
    312             END DO 
    313          END DO 
     273         DO_3D_00_00( 1, jpkm1 ) 
     274            zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     275            z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
     276               &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
     277               &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
     278               &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
     279         END_3D 
    314280         CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    315281         CALL iom_put( "eken", z3d )                 ! kinetic energy 
     
    321287         z3d(1,:, : ) = 0._wp 
    322288         z3d(:,1, : ) = 0._wp 
    323          DO jk = 1, jpkm1 
    324             DO jj = 2, jpj 
    325                DO ji = 2, jpi 
    326                   z3d(ji,jj,jk) = 0.25_wp * ( uu(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)  & 
    327                      &                      + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)  & 
    328                      &                      + vv(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)  & 
    329                      &                      + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)  )  & 
    330                      &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
    331                END DO 
    332             END DO 
    333          END DO 
    334           
     289         DO_3D_00_00( 1, jpkm1 ) 
     290            z3d(ji,jj,jk) = 0.25_wp * ( uu(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)  & 
     291               &                      + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)  & 
     292               &                      + vv(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)  & 
     293               &                      + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)  )  & 
     294               &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
     295         END_3D 
    335296         CALL lbc_lnk( 'diawri', z3d, 'T', 1. ) 
    336297         CALL iom_put( "ke", z3d ) ! kinetic energy 
    337298 
    338299         z2d(:,:)  = 0._wp  
    339          DO jk = 1, jpkm1 
    340             DO jj = 1, jpj 
    341                DO ji = 1, jpi 
    342                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 
    343                END DO 
    344             END DO 
    345          END DO 
     300         DO_3D_11_11( 1, jpkm1 ) 
     301            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 
     302         END_3D 
    346303         CALL iom_put( "ke_zint", z2d )   ! vertically integrated kinetic energy 
    347304 
     
    353310          
    354311         z3d(:,:,jpk) = 0._wp  
    355          DO jk = 1, jpkm1 
    356             DO jj = 1, jpjm1 
    357                DO ji = 1, fs_jpim1   ! vector opt. 
    358                   z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
    359                      &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
    360                END DO 
    361             END DO 
    362          END DO 
     312         DO_3D_00_00( 1, jpkm1 ) 
     313            z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
     314               &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
     315         END_3D 
    363316         CALL lbc_lnk( 'diawri', z3d, 'F', 1. ) 
    364317         CALL iom_put( "relvor", z3d )                  ! relative vorticity 
    365318 
    366          DO jk = 1, jpkm1 
    367             DO jj = 1, jpj 
    368                DO ji = 1, jpi 
    369                   z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)  
    370                END DO 
    371             END DO 
    372          END DO 
     319         DO_3D_11_11( 1, jpkm1 ) 
     320            z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)  
     321         END_3D 
    373322         CALL iom_put( "absvor", z3d )                  ! absolute vorticity 
    374323 
    375          DO jk = 1, jpkm1 
    376             DO jj = 1, jpjm1 
    377                DO ji = 1, fs_jpim1   ! vector opt. 
    378                   ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
    379                      &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    380                   IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
    381                   ELSE                      ;   ze3 = 0._wp 
    382                   ENDIF 
    383                   z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)  
    384                END DO 
    385             END DO 
    386          END DO 
     324         DO_3D_00_00( 1, jpkm1 ) 
     325            ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     326               &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
     327            IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
     328            ELSE                      ;   ze3 = 0._wp 
     329            ENDIF 
     330            z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)  
     331         END_3D 
    387332         CALL lbc_lnk( 'diawri', z3d, 'F', 1. ) 
    388333         CALL iom_put( "potvor", z3d )                  ! potential vorticity 
    389334 
    390335      ENDIF 
    391     
    392336      ! 
    393337      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
     
    404348      IF( iom_use("u_heattr") ) THEN 
    405349         z2d(:,:) = 0._wp  
    406          DO jk = 1, jpkm1 
    407             DO jj = 2, jpjm1 
    408                DO ji = fs_2, fs_jpim1   ! vector opt. 
    409                   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) ) 
    410                END DO 
    411             END DO 
    412          END DO 
     350         DO_3D_00_00( 1, jpkm1 ) 
     351            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) ) 
     352         END_3D 
    413353         CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    414354         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
     
    417357      IF( iom_use("u_salttr") ) THEN 
    418358         z2d(:,:) = 0.e0  
    419          DO jk = 1, jpkm1 
    420             DO jj = 2, jpjm1 
    421                DO ji = fs_2, fs_jpim1   ! vector opt. 
    422                   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) ) 
    423                END DO 
    424             END DO 
    425          END DO 
     359         DO_3D_00_00( 1, jpkm1 ) 
     360            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) ) 
     361         END_3D 
    426362         CALL lbc_lnk( 'diawri', z2d, 'U', -1. ) 
    427363         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
     
    439375      IF( iom_use("v_heattr") ) THEN 
    440376         z2d(:,:) = 0.e0  
    441          DO jk = 1, jpkm1 
    442             DO jj = 2, jpjm1 
    443                DO ji = fs_2, fs_jpim1   ! vector opt. 
    444                   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) ) 
    445                END DO 
    446             END DO 
    447          END DO 
     377         DO_3D_00_00( 1, jpkm1 ) 
     378            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) ) 
     379         END_3D 
    448380         CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    449381         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
     
    452384      IF( iom_use("v_salttr") ) THEN 
    453385         z2d(:,:) = 0._wp  
    454          DO jk = 1, jpkm1 
    455             DO jj = 2, jpjm1 
    456                DO ji = fs_2, fs_jpim1   ! vector opt. 
    457                   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) ) 
    458                END DO 
    459             END DO 
    460          END DO 
     386         DO_3D_00_00( 1, jpkm1 ) 
     387            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) ) 
     388         END_3D 
    461389         CALL lbc_lnk( 'diawri', z2d, 'V', -1. ) 
    462390         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
     
    465393      IF( iom_use("tosmint") ) THEN 
    466394         z2d(:,:) = 0._wp 
    467          DO jk = 1, jpkm1 
    468             DO jj = 2, jpjm1 
    469                DO ji = fs_2, fs_jpim1   ! vector opt. 
    470                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
    471                END DO 
    472             END DO 
    473          END DO 
     395         DO_3D_00_00( 1, jpkm1 ) 
     396            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
     397         END_3D 
    474398         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    475399         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
     
    477401      IF( iom_use("somint") ) THEN 
    478402         z2d(:,:)=0._wp 
    479          DO jk = 1, jpkm1 
    480             DO jj = 2, jpjm1 
    481                DO ji = fs_2, fs_jpim1   ! vector opt. 
    482                   z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
    483                END DO 
    484             END DO 
    485          END DO 
     403         DO_3D_00_00( 1, jpkm1 ) 
     404            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
     405         END_3D 
    486406         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    487407         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
     
    490410      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2) 
    491411      ! 
    492            
     412       
    493413      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
    494414 
     
    506426      INTEGER, DIMENSION(2) :: ierr 
    507427      !!---------------------------------------------------------------------- 
    508       ierr = 0 
    509       ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
    510          &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
    511          &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
     428      IF( nn_write == -1 ) THEN 
     429         dia_wri_alloc = 0 
     430      ELSE     
     431         ierr = 0 
     432         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     & 
     433            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     & 
     434            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) ) 
    512435         ! 
    513       dia_wri_alloc = MAXVAL(ierr) 
    514       CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     436         dia_wri_alloc = MAXVAL(ierr) 
     437         CALL mpp_sum( 'diawri', dia_wri_alloc ) 
     438         ! 
     439      ENDIF 
    515440      ! 
    516441   END FUNCTION dia_wri_alloc 
     442  
     443   INTEGER FUNCTION dia_wri_alloc_abl() 
     444      !!---------------------------------------------------------------------- 
     445     ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl) 
     446      CALL mpp_sum( 'diawri', dia_wri_alloc_abl ) 
     447      ! 
     448   END FUNCTION dia_wri_alloc_abl 
    517449 
    518450    
     
    538470      INTEGER  ::   ierr                                     ! error code return from allocation 
    539471      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers 
     472      INTEGER  ::   ipka                                     ! ABL 
    540473      INTEGER  ::   jn, ierror                               ! local integers 
    541474      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars 
     
    543476      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace 
    544477      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace 
     478      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace 
    545479      !!---------------------------------------------------------------------- 
    546480      ! 
     
    576510      ijmi = 1      ;      ijma = jpj 
    577511      ipk = jpk 
     512      IF(ln_abl) ipka = jpkam1 
    578513 
    579514      ! define time axis 
     
    678613            &          "m", ipk, gdepw_1d, nz_W, "down" ) 
    679614 
     615         IF( ln_abl ) THEN  
     616         ! Define the ABL grid FILE ( nid_A ) 
     617            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' ) 
     618            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename 
     619            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
     620               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
     621               &          nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set ) 
     622            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept 
     623               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" ) 
     624            !                                                            ! Index of ocean points 
     625         ALLOCATE( zw3d_abl(jpi,jpj,ipka) )  
     626         zw3d_abl(:,:,:) = 1._wp  
     627         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume 
     628            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface 
     629         DEALLOCATE(zw3d_abl) 
     630         ENDIF 
    680631 
    681632         ! Declare all the output fields as NETCDF variables 
     
    727678         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm 
    728679            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    729 ! 
     680         ! 
     681         IF( ln_abl ) THEN 
     682            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl 
     683               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     684            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl 
     685               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     686            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl 
     687               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
     688            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl 
     689               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     690            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl 
     691               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     692            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl 
     693               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     694            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl 
     695               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )  
     696            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh 
     697               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                  
     698#if defined key_si3 
     699            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i 
     700               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout ) 
     701#endif 
     702            CALL histend( nid_A, snc4chunks=snc4set ) 
     703         ENDIF 
     704         ! 
    730705         IF( ln_icebergs ) THEN 
    731706            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , & 
     
    885860      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction    
    886861      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed    
    887 ! 
     862      ! 
     863      IF( ln_abl ) THEN  
     864         ALLOCATE( zw3d_abl(jpi,jpj,jpka) ) 
     865         IF( ln_mskland )   THEN  
     866            DO jk=1,jpka 
     867               zw3d_abl(:,:,jk) = tmask(:,:,1) 
     868            END DO        
     869         ELSE 
     870            zw3d_abl(:,:,:) = 1._wp      
     871         ENDIF        
     872         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh  
     873         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl 
     874         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl 
     875         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl 
     876         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl        
     877         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl 
     878         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl 
     879         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl  
     880#if defined key_si3 
     881         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i 
     882#endif 
     883         DEALLOCATE(zw3d_abl) 
     884      ENDIF 
     885      ! 
    888886      IF( ln_icebergs ) THEN 
    889887         ! 
     
    931929      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    932930 
    933       CALL histwrite( nid_W, "vovecrtz", it, ww             , ndim_T, ndex_T )    ! vert. current 
     931      IF( ln_zad_Aimp ) THEN 
     932         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current 
     933      ELSE 
     934         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current 
     935      ENDIF 
    934936      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    935937      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     
    951953         CALL histclo( nid_V ) 
    952954         CALL histclo( nid_W ) 
     955         IF(ln_abl) CALL histclo( nid_A ) 
    953956      ENDIF 
    954957      ! 
     
    974977      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created 
    975978      !! 
    976       INTEGER :: inum 
     979      INTEGER :: inum, jk 
    977980      !!---------------------------------------------------------------------- 
    978981      !  
     
    981984      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created ' 
    982985      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc' 
    983  
    984 #if defined key_si3 
    985      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) 
    986 #else 
    987      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
    988 #endif 
    989  
     986      ! 
     987      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. ) 
     988      ! 
    990989      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature 
    991990      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity 
     
    993992      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity 
    994993      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity 
    995       CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww                )    ! now k-velocity 
     994      IF( ln_zad_Aimp ) THEN 
     995         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity 
     996      ELSE 
     997         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity 
     998      ENDIF 
     999      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity 
     1000      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height 
     1001      ! 
     1002      IF ( ln_isf ) THEN 
     1003         IF (ln_isfcav_mlt) THEN 
     1004            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity 
     1005            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity 
     1006            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity 
     1007            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity 
     1008            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity 
     1009            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 
     1010         END IF 
     1011         IF (ln_isfpar_mlt) THEN 
     1012            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity 
     1013            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity 
     1014            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity 
     1015            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity 
     1016            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity 
     1017            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity 
     1018            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 
     1019         END IF 
     1020      END IF 
     1021      ! 
    9961022      IF( ALLOCATED(ahtu) ) THEN 
    9971023         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
     
    10171043         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity 
    10181044      ENDIF 
    1019   
     1045      IF ( ln_abl ) THEN 
     1046         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind 
     1047         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind 
     1048         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature 
     1049         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity 
     1050      ENDIF 
     1051      ! 
     1052      CALL iom_close( inum ) 
     1053      !  
    10201054#if defined key_si3 
    10211055      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid 
     1056         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' ) 
    10221057         CALL ice_wri_state( inum ) 
     1058         CALL iom_close( inum ) 
    10231059      ENDIF 
    10241060#endif 
    1025       ! 
    1026       CALL iom_close( inum ) 
    1027       !  
     1061 
    10281062   END SUBROUTINE dia_wri_state 
    10291063 
Note: See TracChangeset for help on using the changeset viewer.