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 8055 for branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90 – NEMO

Ignore:
Timestamp:
2017-05-20T13:49:38+02:00 (7 years ago)
Author:
gm
Message:

wrk_OMP: 2nd step: add OMP processor distribution in ZDF package

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90

    r7990 r8055  
    3535   PUBLIC   zdf_iwm         ! called in step module  
    3636   PUBLIC   zdf_iwm_init    ! called in nemogcm module  
    37    PUBLIC   zdf_iwm_alloc   ! called in nemogcm module 
    3837 
    3938   !                       !!* Namelist  namzdf_iwm : internal wave-driven mixing * 
     
    5655 
    5756   !! * Substitutions 
    58 #  include "vectopt_loop_substitute.h90" 
     57#  include "domain_substitute.h90"    
    5958   !!---------------------------------------------------------------------- 
    6059   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     
    7877 
    7978 
    80    SUBROUTINE zdf_iwm( kt ) 
     79   SUBROUTINE zdf_iwm( ARG_2D, kt, p_avm, p_avt, p_avs ) 
    8180      !!---------------------------------------------------------------------- 
    8281      !!                  ***  ROUTINE zdf_iwm  *** 
     
    127126      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
    128127      !!---------------------------------------------------------------------- 
    129       INTEGER, INTENT(in) ::   kt   ! ocean time-step  
     128      INTEGER                    , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices 
     129      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     130      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points) 
     131      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points) 
    130132      ! 
    131133      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    132       REAL(wp) ::   ztpc         ! scalar workspace 
    133       REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure 
    134       REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth 
    135       REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom 
    136       REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution 
    137       REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid) 
    138       REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid) 
    139       REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter 
     134      REAL(wp) ::   zztemp       ! scalar workspace 
     135      REAL(wp), DIMENSION(WRK_2D) ::  zfact     ! Used for vertical structure 
     136      REAL(wp), DIMENSION(WRK_2D) ::  zhdep     ! Ocean depth 
     137      REAL(wp), DIMENSION(WRK_3D) ::  zwkb      ! WKB-stretched height above bottom 
     138      REAL(wp), DIMENSION(WRK_3D) ::  zweight   ! Weight for high mode vertical distribution 
     139      REAL(wp), DIMENSION(WRK_3D) ::  znu_t     ! Molecular kinematic viscosity (T grid) 
     140      REAL(wp), DIMENSION(WRK_3D) ::  znu_w     ! Molecular kinematic viscosity (W grid) 
     141      REAL(wp), DIMENSION(WRK_3D) ::  zReb      ! Turbulence intensity parameter 
    140142      !!---------------------------------------------------------------------- 
    141143      ! 
    142144      IF( nn_timing == 1 )   CALL timing_start('zdf_iwm') 
    143145      ! 
    144       CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
    145       CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    146  
    147146      !                          ! ----------------------------- ! 
    148147      !                          !  Internal wave-driven mixing  !  (compute zav_wave) 
     
    151150      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    152151      !                                                 using an exponential decay from the seafloor. 
    153       DO jj = 1, jpj                ! part independent of the level 
    154          DO ji = 1, jpi 
     152      DO jj = k_Jstr, k_Jend 
     153         DO ji = k_Istr, k_Iend     ! part independent of the level 
    155154            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    156155            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     
    158157         END DO 
    159158      END DO 
    160  
    161159      DO jk = 2, jpkm1              ! complete with the level-dependent part 
    162          emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      & 
    163             &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   & 
     160         DO jj = k_Jstr, k_Jend 
     161            DO ji = k_Istr, k_Iend 
     162               emix_iwm(ji,jj,jk) = zfact(ji,jj) * wmask(ji,jj,jk)                              &  
     163                  &   * (  EXP( ( gde3w_n(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )      & 
     164                  &      - EXP( ( gde3w_n(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )  )   & 
     165                  &   / ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) 
     166 
    164167!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point 
    165             &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
     168!!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all 
     169 
     170            END DO 
     171         END DO 
    166172      END DO 
    167173 
    168174      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
    169175      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
    170       ! 
     176      !                                          ! (NB: N2 is masked, so no use of wmask here) 
    171177      SELECT CASE ( nn_zpyc ) 
    172178      ! 
    173179      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    174180         ! 
    175          zfact(:,:) = 0._wp 
    176          DO jk = 2, jpkm1              ! part independent of the level 
    177             zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    178          END DO 
    179          ! 
    180          DO jj = 1, jpj 
    181             DO ji = 1, jpi 
    182                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    183             END DO 
    184          END DO 
     181         zfact(WRK_2D) = 0._wp         ! part independent of the level 
     182         DO jk = 2, jpkm1 
     183            zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  ) 
     184         END DO 
     185         ! 
     186         WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 
     187!         DO jj = k_Jstr, k_Jend 
     188!            DO ji = k_Istr, k_Iend 
     189!               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     190!            END DO 
     191!         END DO 
     192         !                             ! complete with the level-dependent part 
     193         DO jk = 2, jpkm1 
     194            emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zfact(WRK_2D) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  ) 
     195         END DO 
     196         ! 
     197      CASE ( 2 )               ! Dissipation scales as N^2 
     198         ! 
     199         zfact(WRK_2D) = 0._wp 
     200         DO jk = 2, jpkm1              ! part independent of the level  
     201            zfact(WRK_2D) = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * MAX( 0._wp, rn2(WRK_2D,jk) ) 
     202         END DO 
     203         ! 
     204         WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = epyc_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 
     205!         DO jj = k_Jstr, k_Jend 
     206!            DO ji = k_Istr, k_Iend 
     207!               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     208!            END DO 
     209!         END DO 
    185210         ! 
    186211         DO jk = 2, jpkm1              ! complete with the level-dependent part 
    187             emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    188          END DO 
    189          ! 
    190       CASE ( 2 )               ! Dissipation scales as N^2 
    191          ! 
    192          zfact(:,:) = 0._wp 
    193          DO jk = 2, jpkm1              ! part independent of the level 
    194             zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    195          END DO 
    196          ! 
    197          DO jj= 1, jpj 
    198             DO ji = 1, jpi 
    199                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    200             END DO 
    201          END DO 
    202          ! 
    203          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    204             emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     212            emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zfact(WRK_2D) * MAX( 0._wp, rn2(WRK_2D,jk) ) 
    205213         END DO 
    206214         ! 
     
    210218      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    211219      ! 
    212       zwkb(:,:,:) = 0._wp 
    213       zfact(:,:) = 0._wp 
     220      zwkb (WRK_3D) = 0._wp 
     221      zfact(WRK_2D) = 0._wp 
    214222      DO jk = 2, jpkm1 
    215          zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    216          zwkb(:,:,jk) = zfact(:,:) 
     223!!gm  this is potentially dangerous  
     224!         zfact(WRK_2D)    = zfact(WRK_2D) + e3w_n(WRK_2D,jk) * SQRT(  MAX( 0._wp, rn2(WRK_2D,jk) )  ) 
     225!         zwkb (WRK_2D,jk) = zfact(WRK_2D) 
     226         DO jj = k_Jstr, k_Jend 
     227            DO ji = k_Istr, k_Iend 
     228               zztemp = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) 
     229               zfact(ji,jj)    = zztemp 
     230               zwkb (ji,jj,jk) = zztemp 
     231            END DO 
     232         END DO          
    217233      END DO 
    218234      ! 
    219235      DO jk = 2, jpkm1 
    220          DO jj = 1, jpj 
    221             DO ji = 1, jpi 
     236         DO jj = k_Jstr, k_Jend 
     237            DO ji = k_Istr, k_Iend 
    222238               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
    223                   &                                     * tmask(ji,jj,jk) / zfact(ji,jj) 
     239                  &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
    224240            END DO 
    225241         END DO 
    226242      END DO 
    227       zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
    228       ! 
    229       zweight(:,:,:) = 0._wp 
     243      zwkb(WRK_2D,1) = zhdep(WRK_2D) * wmask(WRK_2D,1) 
     244      ! 
     245      zweight(WRK_3D) = 0._wp 
    230246      DO jk = 2, jpkm1 
    231          zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    & 
    232             &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  ) 
    233       END DO 
    234       ! 
    235       zfact(:,:) = 0._wp 
     247         zweight(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * hbot_iwm(:,:)     & 
     248            &               * (  EXP( -zwkb(WRK_2D,jk  ) / hbot_iwm(WRK_2D) )  & 
     249            &                  - EXP( -zwkb(WRK_2D,jk-1) / hbot_iwm(WRK_2D) )  ) 
     250      END DO 
     251      ! 
     252      zfact(WRK_2D) = 0._wp 
    236253      DO jk = 2, jpkm1              ! part independent of the level 
    237          zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    238       END DO 
    239       ! 
    240       DO jj = 1, jpj 
    241          DO ji = 1, jpi 
    242             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    243          END DO 
    244       END DO 
     254         zfact(WRK_2D) = zfact(WRK_2D) + zweight(WRK_2D,jk) 
     255      END DO 
     256      ! 
     257      WHERE( zfact(WRK_2D) /= 0 )   zfact(WRK_2D) = ebot_iwm(WRK_2D) / ( rau0 * zfact(WRK_2D) ) 
     258!      DO jj = k_Jstr, k_Jend 
     259!         DO ji = k_Istr, k_Iend 
     260!            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     261!         END DO 
     262!      END DO 
    245263      ! 
    246264      DO jk = 2, jpkm1              ! complete with the level-dependent part 
    247          emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    248             &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    249       END DO 
    250       ! 
     265         emix_iwm(WRK_2D,jk) = emix_iwm(WRK_2D,jk) + zweight(WRK_2D,jk) * zfact(WRK_2D) * wmask(WRK_2D,jk)   & 
     266            &                                / ( gde3w_n(WRK_2D,jk) - gde3w_n(WRK_2D,jk-1) ) 
     267!!gm  use of e3t_n just above? 
     268      END DO 
     269      ! 
     270!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    251271      ! Calculate molecular kinematic viscosity 
    252       znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
    253          &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
     272      znu_t(WRK_3D) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(WRK_3D,jp_tem)             & 
     273         &                        + 0.00694_wp * tsn(WRK_3D,jp_tem) * tsn(WRK_3D,jp_tem)   & 
     274         &                        + 0.02305_wp * tsn(WRK_3D,jp_sal)  ) * tmask(WRK_3D) * r1_rau0 
    254275      DO jk = 2, jpkm1 
    255          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    256       END DO 
     276         znu_w(WRK_2D,jk) = 0.5_wp * ( znu_t(WRK_2D,jk-1) + znu_t(WRK_2D,jk) ) * wmask(WRK_2D,jk) 
     277      END DO 
     278!!gm end 
    257279      ! 
    258280      ! Calculate turbulence intensity parameter Reb 
    259281      DO jk = 2, jpkm1 
    260          zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
     282         zReb(WRK_2D,jk) = emix_iwm(WRK_2D,jk) / MAX( 1.e-20_wp, znu_w(WRK_2D,jk) * rn2(WRK_2D,jk) ) 
    261283      END DO 
    262284      ! 
    263285      ! Define internal wave-induced diffusivity 
    264286      DO jk = 2, jpkm1 
    265          zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     287         zav_wave(WRK_2D,jk) = znu_w(WRK_2D,jk) * zReb(WRK_2D,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    266288      END DO 
    267289      ! 
    268290      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    269291         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
    270             DO jj = 1, jpj 
    271                DO ji = 1, jpi 
     292            DO jj = k_Jstr, k_Jend 
     293               DO ji = k_Istr, k_Iend 
    272294                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    273295                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     
    281303      ! 
    282304      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    283          zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
     305         zav_wave(WRK_2D,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(WRK_2D,jk) ), 1.e-2_wp  ) * wmask(WRK_2D,jk) 
    284306      END DO 
    285307      ! 
    286308      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    287          ztpc = 0._wp 
     309         zztemp = 0._wp 
    288310!!gm used of glosum 3D.... 
    289311         DO jk = 2, jpkm1 
    290             DO jj = 1, jpj 
    291                DO ji = 1, jpi 
    292                   ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
     312            DO jj = k_Jstr, k_Jend 
     313               DO ji = k_Istr, k_Iend 
     314                  zztemp = zztemp + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
    293315                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    294316               END DO 
    295317            END DO 
    296318         END DO 
    297          IF( lk_mpp )   CALL mpp_sum( ztpc ) 
    298          ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     319         IF( lk_mpp )   CALL mpp_sum( zztemp ) 
     320         zztemp = rau0 * zztemp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
    299321         ! 
    300322         IF(lwp) THEN 
     
    303325            WRITE(numout,*) '~~~~~~~ ' 
    304326            WRITE(numout,*) 
    305             WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW' 
     327            WRITE(numout,*) '      Total power consumption by av_wave =  ', zztemp * 1.e-12_wp, 'TW' 
    306328         ENDIF 
    307329      ENDIF 
     
    313335      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
    314336         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb 
    315             DO jj = 1, jpj 
    316                DO ji = 1, jpi 
     337            DO jj = k_Jstr, k_Jend 
     338               DO ji = k_Istr, k_Iend 
    317339                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  & 
    318340                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   & 
     
    323345         CALL iom_put( "av_ratio", zav_ratio ) 
    324346         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
    325             avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
    326             avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
    327             avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     347            p_avs(WRK_2D,jk) = p_avs(WRK_2D,jk) + zav_wave(WRK_2D,jk) * zav_ratio(WRK_2D,jk) 
     348            p_avt(WRK_2D,jk) = p_avt(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
     349            p_avm(WRK_2D,jk) = p_avm(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
    328350         END DO 
    329351         ! 
    330352      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    331353         DO jk = 2, jpkm1 
    332             avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) 
    333             avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
    334             avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     354            p_avs(WRK_2D,jk) = p_avs(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
     355            p_avt(WRK_2D,jk) = p_avt(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
     356            p_avm(WRK_2D,jk) = p_avm(WRK_2D,jk) + zav_wave(WRK_2D,jk) 
    335357         END DO 
    336358      ENDIF 
     
    341363                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm) 
    342364      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    343          bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    344          pcmap_iwm(:,:) = 0._wp 
    345365         DO jk = 2, jpkm1 
    346             pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk) 
    347          END DO 
    348          pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:) 
    349          CALL iom_put( "bflx_iwm", bflx_iwm ) 
     366            bflx_iwm(WRK_2D,jk) = MAX( 0._wp, rn2(WRK_2D,jk) ) * zav_wave(WRK_2D,jk) 
     367         END DO 
     368         pcmap_iwm(WRK_2D) = 0._wp 
     369         DO jk = 2, jpkm1 
     370            pcmap_iwm(WRK_2D) = pcmap_iwm(WRK_2D) + e3w_n(WRK_2D,jk) * bflx_iwm(WRK_2D,jk) 
     371         END DO 
     372         pcmap_iwm(WRK_2D) = rau0 * pcmap_iwm(WRK_2D) 
     373         CALL iom_put( "bflx_iwm" , bflx_iwm  ) 
    350374         CALL iom_put( "pcmap_iwm", pcmap_iwm ) 
    351375      ENDIF 
     
    353377      CALL iom_put( "emix_iwm", emix_iwm ) 
    354378       
    355       CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
    356       CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
    357  
    358379      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
    359380      ! 
Note: See TracChangeset for help on using the changeset viewer.