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 15517 – NEMO

Changeset 15517


Ignore:
Timestamp:
2021-11-16T15:48:25+01:00 (6 months ago)
Author:
cdllod
Message:

zdfiwm optimization by Clement

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r15388_updated_zdfiwm/src/OCE/ZDF/zdfiwm.F90

    r15466 r15517  
    4242   LOGICAL ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F) 
    4343 
    44    REAL(wp)::  r1_6 = 1._wp / 6._wp 
    45    REAL(wp)::  rnu  = 1.4e-6_wp   ! molecular kinematic viscosity 
     44   REAL(wp)::  z1_6 = 1._wp / 6._wp 
     45   REAL(wp)::  znu  = 1.4e-6_wp   ! molecular kinematic viscosity 
    4646 
    4747   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ebot_iwm   ! bottom-intensified dissipation above abyssal hills (W/m2) 
     
    135135      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    136136      REAL(wp), SAVE :: zztmp 
    137       REAL(wp)       :: ztmp1, ztmp2        ! scalar workspace 
     137      ! 
    138138      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zfact       ! Used for vertical structure 
    139139      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zReb        ! Turbulence intensity parameter 
     
    145145      !!---------------------------------------------------------------------- 
    146146      ! 
    147       !                       !* Set to zero the 1st and last vertical levels of appropriate variables 
    148       IF( iom_use("emix_iwm") ) THEN 
    149          zemx_iwm(:,:,:) = 0._wp 
    150       ENDIF 
    151       IF( iom_use("av_ratio") ) THEN 
    152          zav_ratio(:,:,:) = 0._wp 
    153       ENDIF 
    154       IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN 
    155          zav_wave(:,:,:) = 0._wp 
    156       ENDIF 
     147      !                       !* Initialize appropriately certain variables 
     148      zav_ratio(:,:,:) = 1._wp * wmask(:,:,:)  ! important to set it to 1 here  
     149      IF( iom_use("emix_iwm") )                         zemx_iwm (:,:,:) = 0._wp 
     150      IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl )   zav_wave (:,:,:) = 0._wp 
    157151      ! 
    158152      !                       ! ----------------------------- ! 
     
    163157      !                       !* ocean depth using an exponential decay from the seafloor. 
    164158      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                ! part independent of the level 
    165          zfact(ji,jj) = rho0 * (  1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) )  ) 
    166          IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 
     159         IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ecri_iwm(ji,jj) * r1_rho0 / ( 1._wp - EXP( -ht(ji,jj) * hcri_iwm(ji,jj) ) ) 
     160         ELSE                          ; zfact(ji,jj) = 0._wp 
     161         ENDIF 
    167162      END_2D 
    168163 
    169164      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )      ! complete with the level-dependent part 
    170          IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
    171             zemx_iwm(ji,jj,jk) = 0._wp 
    172          ELSE 
    173             zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gdept(ji,jj,jk  ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) )     & 
    174                  &                               - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) ) )   & 
    175                  &                            / e3w(ji,jj,jk,Kmm) 
    176          ENDIF 
     165         zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gdept(ji,jj,jk  ,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) )   & 
     166            &                                 - EXP( ( gdept(ji,jj,jk-1,Kmm) - ht(ji,jj) ) * hcri_iwm(ji,jj) )   & 
     167            &                                ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm) 
    177168      END_3D 
    178169 
    179170                               !* 'bot' component: distribute energy over the time-varying 
    180171                               !* ocean depth using an algebraic decay above the seafloor. 
    181       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    182          zfact(ji,jj) = 0._wp 
    183       END_2D 
    184172      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )               ! part independent of the level 
    185          IF( ht(ji,jj) /= 0._wp )  zfact(ji,jj) = ebot_iwm(ji,jj) * (  1._wp +  hbot_iwm(ji,jj) / ht(ji,jj)  ) / rho0 
     173         IF( ht(ji,jj) /= 0._wp ) THEN ; zfact(ji,jj) = ebot_iwm(ji,jj) * ( 1._wp +  hbot_iwm(ji,jj) / ht(ji,jj)  ) * r1_rho0 
     174         ELSE                          ; zfact(ji,jj) = 0._wp 
     175         ENDIF 
    186176      END_2D 
    187177 
    188178      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )     ! complete with the level-dependent part 
    189          zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) +                                                                                          &  
    190             &         zfact(ji,jj) * (   1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk  ,Kmm) ) / hbot_iwm(ji,jj) )                        & 
    191             &                         -  1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) )   ) * wmask(ji,jj,jk)  & 
    192             &                      / e3w(ji,jj,jk,Kmm) 
     179         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) +                                                                           &  
     180            &                 zfact(ji,jj) * (  1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk  ,Kmm) ) / hbot_iwm(ji,jj) )  & 
     181            &                                 - 1._wp / ( 1._wp + ( ht(ji,jj) - gdept(ji,jj,jk-1,Kmm) ) / hbot_iwm(ji,jj) )  & 
     182            &                                ) * wmask(ji,jj,jk) / e3w(ji,jj,jk,Kmm) 
    193183      END_3D 
    194184 
     
    203193      ! 
    204194      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    205          IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ensq_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
     195         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ensq_iwm(ji,jj) * r1_rho0 / zfact(ji,jj) 
    206196      END_2D 
    207197      ! 
     
    216206      END_2D 
    217207      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )     ! part independent of the level 
    218          zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) ) ) 
     208         zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) 
    219209      END_3D 
    220210      ! 
    221211      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    222          IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = esho_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
     212         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = esho_iwm(ji,jj) * r1_rho0 / zfact(ji,jj) 
    223213      END_2D 
    224214      ! 
    225215      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )     ! complete with the level-dependent part 
    226          zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) ) ) 
     216         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) 
    227217      END_3D 
    228218 
    229219      ! Calculate turbulence intensity parameter Reb 
    230220      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    231          zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, rnu * rn2(ji,jj,jk) ) 
     221         zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu * rn2(ji,jj,jk) ) 
    232222      END_3D 
    233223      ! 
    234224      ! Define internal wave-induced diffusivity 
    235225      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    236          zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * r1_6 * rnu  ! This corresponds to a constant mixing efficiency of 1/6 
     226         zav_wave(ji,jj,jk) = zReb(ji,jj,jk) * z1_6 * znu  ! This corresponds to a constant mixing efficiency of 1/6 
    237227      END_3D 
    238228      ! 
     
    240230         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224) regimes 
    241231            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    242                zav_wave(ji,jj,jk) = 3.6515_wp * rnu * SQRT( zReb(ji,jj,jk) ) 
     232               zav_wave(ji,jj,jk) = 3.6515_wp * znu * SQRT( zReb(ji,jj,jk) ) 
    243233            ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 
    244                zav_wave(ji,jj,jk) = 0.052125_wp * rnu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     234               zav_wave(ji,jj,jk) = 0.052125_wp * znu * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
    245235            ENDIF 
    246236         END_3D 
     
    248238      ! 
    249239      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )    ! Bound diffusivity by molecular value and 100 cm2/s 
    250          zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk) 
    251       END_3D 
     240         zav_wave(ji,jj,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp ) * wmask(ji,jj,jk) 
     241      END_3D 
     242      ! 
     243      !                          ! ----------------------- ! 
     244      !                          !   Update  mixing coefs  !                           
     245      !                          ! ----------------------- ! 
     246      ! 
     247      IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature 
     248         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb (else it is set to 1) 
     249            zav_ratio(ji,jj,jk) = ( 0.505_wp + & 
     250               &                    0.495_wp * TANH( 0.92_wp * ( LOG10( MAX( 1.e-20, zReb(ji,jj,jk) * 5._wp * z1_6 ) ) - 0.60_wp ) ) & 
     251               &                  ) * wmask(ji,jj,jk) 
     252         END_3D 
     253      ENDIF 
     254      CALL iom_put( "av_ratio", zav_ratio ) 
     255      ! 
     256      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing 
     257         p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 
     258         p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     259         p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     260      END_3D 
     261      !                             !* output internal wave-driven mixing coefficient 
     262      CALL iom_put( "av_wave", zav_wave ) 
     263                                    !* output useful diagnostics: Kz*N^2 ,  
     264                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
     265      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
     266         ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) ) 
     267         z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp   ! Initialisation for iom_put 
     268         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     269            z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) 
     270            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk) 
     271         END_3D 
     272         CALL iom_put(  "bflx_iwm", z3d ) 
     273         CALL iom_put( "pcmap_iwm", z2d ) 
     274         DEALLOCATE( z2d , z3d ) 
     275      ENDIF 
     276      CALL iom_put( "emix_iwm", zemx_iwm ) 
     277 
    252278      ! 
    253279      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
     
    260286         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile 
    261287            CALL mpp_sum( 'zdfiwm', zztmp ) 
    262             zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
     288            zztmp = rho0 * zztmp ! Global integral of rho0 * Kz * N^2 = power contributing to mixing 
    263289            ! 
    264290            IF(lwp) THEN 
     
    272298      ENDIF 
    273299 
    274       !                          ! ----------------------- ! 
    275       !                          !   Update  mixing coefs  !                           
    276       !                          ! ----------------------- ! 
    277       ! 
    278       IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature 
    279          ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 
    280          DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb 
    281             ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 
    282             IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 
    283                zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) ) 
    284             ELSE 
    285                zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk) 
    286             ENDIF 
    287          END_3D 
    288          CALL iom_put( "av_ratio", zav_ratio ) 
    289          DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) !* update momentum & tracer diffusivity with wave-driven mixing 
    290             p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 
    291             p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
    292             p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
    293          END_3D 
    294          ! 
    295       ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing 
    296          DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    297             p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) 
    298             p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
    299             p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
    300          END_3D 
    301       ENDIF 
    302       !                             !* output internal wave-driven mixing coefficient 
    303       CALL iom_put( "av_wave", zav_wave ) 
    304                                     !* output useful diagnostics: Kz*N^2 ,  
    305                                     !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
    306       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    307          ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) ) 
    308          ! Initialisation for iom_put 
    309          z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp 
    310  
    311          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    312             z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) 
    313             z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk) 
    314          END_3D 
    315          DO_2D( 0, 0, 0, 0 ) 
    316             z2d(ji,jj) = rho0 * z2d(ji,jj) 
    317          END_2D 
    318          CALL iom_put(  "bflx_iwm", z3d ) 
    319          CALL iom_put( "pcmap_iwm", z2d ) 
    320          DEALLOCATE( z2d , z3d ) 
    321       ENDIF 
    322       CALL iom_put( "emix_iwm", zemx_iwm ) 
    323  
    324300      IF(sn_cfctl%l_prtctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk) 
    325301      ! 
     
    356332      INTEGER  ::   inum               ! local integer 
    357333      INTEGER  ::   ios 
    358       REAL(wp) ::   zbot, zcri, znsq, zsho   ! local scalars 
    359334      ! 
    360335      CHARACTER(len=256)            ::   cn_dir                 ! Root directory for location of ssr files 
     
    372347      TYPE(FLD  ), DIMENSION(jpiwm) ::   sf_iwm                         ! structure of input fields (file informations, fields read) 
    373348      ! 
     349      REAL(wp), DIMENSION(jpi,jpj,4) ::   ztmp 
     350      REAL(wp), DIMENSION(4)         ::   zdia 
     351      ! 
    374352      NAMELIST/namzdf_iwm/ ln_mevar, ln_tsdiff, & 
    375353          &                cn_dir, sn_mpb, sn_mpc, sn_mpn, sn_mps, sn_dsb, sn_dsc 
     
    395373      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should  
    396374      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 
    397       avmb(:) = rnu              ! molecular value 
     375      avmb(:) = znu              ! molecular value 
    398376      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_iwm)     
    399       avtb_2d(:,:) = 1._wp     ! uniform  
     377      avtb_2d(:,:) = 1._wp       ! uniform  
    400378      IF(lwp) THEN                  ! Control print 
    401379         WRITE(numout,*) 
     
    422400      sf_iwm(jp_mpb)%fnow(:,:,1) = 1.e-10_wp 
    423401      sf_iwm(jp_mpc)%fnow(:,:,1) = 1.e-10_wp 
    424       sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e-6_wp 
     402      sf_iwm(jp_mpn)%fnow(:,:,1) = 1.e-5_wp 
    425403      sf_iwm(jp_mps)%fnow(:,:,1) = 1.e-10_wp 
    426404      sf_iwm(jp_dsb)%fnow(:,:,1) = 100._wp 
     
    439417      hcri_iwm(:,:) = 1._wp / hcri_iwm(:,:) ! only the inverse height is used, hence calculated here once for all 
    440418 
    441       zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) ) 
    442       zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) ) 
    443       znsq = glob_sum( 'zdfiwm', e1e2t(:,:) * ensq_iwm(:,:) ) 
    444       zsho = glob_sum( 'zdfiwm', e1e2t(:,:) * esho_iwm(:,:) ) 
     419      ! diags 
     420      ztmp(:,:,1) = e1e2t(:,:) * ebot_iwm(:,:) 
     421      ztmp(:,:,2) = e1e2t(:,:) * ecri_iwm(:,:) 
     422      ztmp(:,:,3) = e1e2t(:,:) * ensq_iwm(:,:) 
     423      ztmp(:,:,4) = e1e2t(:,:) * esho_iwm(:,:) 
     424 
     425      zdia(1:4) = glob_sum_vec( 'zdfiwm', ztmp(:,:,1:4) ) 
    445426 
    446427      IF(lwp) THEN 
    447          WRITE(numout,*) '      Dissipation above abyssal hills:        ', zbot * 1.e-12_wp, 'TW' 
    448          WRITE(numout,*) '      Dissipation along topographic slopes:   ', zcri * 1.e-12_wp, 'TW' 
    449          WRITE(numout,*) '      Dissipation scaling with N^2:           ', znsq * 1.e-12_wp, 'TW' 
    450          WRITE(numout,*) '      Dissipation due to shoaling:            ', zsho * 1.e-12_wp, 'TW' 
     428         WRITE(numout,*) '      Dissipation above abyssal hills:        ', zdia(1) * 1.e-12_wp, 'TW' 
     429         WRITE(numout,*) '      Dissipation along topographic slopes:   ', zdia(2) * 1.e-12_wp, 'TW' 
     430         WRITE(numout,*) '      Dissipation scaling with N^2:           ', zdia(3) * 1.e-12_wp, 'TW' 
     431         WRITE(numout,*) '      Dissipation due to shoaling:            ', zdia(4) * 1.e-12_wp, 'TW' 
    451432      ENDIF 
    452433      ! 
Note: See TracChangeset for help on using the changeset viewer.