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 13463 for NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/ZDF/zdfiwm.F90 – NEMO

Ignore:
Timestamp:
2020-09-14T17:40:34+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2195:update to trunk 13461

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

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
         8 
         9# SETTE 
         10^/utils/CI/sette@13382        sette 

  • NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/ZDF/zdfiwm.F90

    r10425 r13463  
    2323   USE phycst         ! physical constants 
    2424   ! 
     25   USE fldread        ! field read 
    2526   USE prtctl         ! Print control 
    2627   USE in_out_manager ! I/O manager 
    … …  
    4950 
    5051   !! * Substitutions 
    51 #  include "vectopt_loop_substitute.h90" 
     52#  include "do_loop_substitute.h90" 
     53#  include "domzgr_substitute.h90" 
    5254   !!---------------------------------------------------------------------- 
    5355   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    6971 
    7072 
    71    SUBROUTINE zdf_iwm( kt, p_avm, p_avt, p_avs ) 
     73   SUBROUTINE zdf_iwm( kt, Kmm, p_avm, p_avt, p_avs ) 
    7274      !!---------------------------------------------------------------------- 
    7375      !!                  ***  ROUTINE zdf_iwm  *** 
    … …  
    8789      !!              This is divided into three components: 
    8890      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
    89       !!                     zemx_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm ) 
     91      !!                     zemx_iwm(z) = ( ecri_iwm / rho0 ) * EXP( -(H-z)/hcri_iwm ) 
    9092      !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm 
    9193      !!              where hcri_iwm is the characteristic length scale of the bottom  
    9294      !!              intensification, ecri_iwm a map of available power, and H the ocean depth. 
    9395      !!                 2. Pycnocline-intensified low-mode dissipation 
    94       !!                     zemx_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
    95       !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
     96      !!                     zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     97      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w[z) ) 
    9698      !!              where epyc_iwm is a map of available power, and nn_zpyc 
    9799      !!              is the chosen stratification-dependence of the internal wave 
    98100      !!              energy dissipation. 
    99101      !!                 3. WKB-height dependent high mode dissipation 
    100       !!                     zemx_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
    101       !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 
     102      !!                     zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
     103      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w[z) ) 
    102104      !!              where hbot_iwm is the characteristic length scale of the WKB bottom  
    103105      !!              intensification, ebot_iwm is a map of available power, and z_wkb is the 
    104106      !!              WKB-stretched height above bottom defined as 
    105       !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) 
    106       !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    ) 
     107      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w[z'>=z) ) 
     108      !!                                 / SUM( sqrt(rn2(z'))    * e3w[z')    ) 
    107109      !! 
    108110      !!              - update the model vertical eddy viscosity and diffusivity:  
    … …  
    118120      !!---------------------------------------------------------------------- 
    119121      INTEGER                    , INTENT(in   ) ::   kt             ! ocean time step 
     122      INTEGER                    , INTENT(in   ) ::   Kmm            ! time level index 
    120123      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avm          ! momentum Kz (w-points) 
    121124      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   p_avt, p_avs   ! tracer   Kz (w-points) 
    … …  
    137140      !!---------------------------------------------------------------------- 
    138141      ! 
    139       !                       !* Set to zero the 1st and last vertical levels of appropriate variables 
    140       zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp 
    141       zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp 
    142       zav_wave (:,:,1) = 0._wp   ;   zav_wave (:,:,jpk) = 0._wp 
     142      !                        
     143      ! Set to zero the 1st and last vertical levels of appropriate variables 
     144      IF( iom_use("emix_iwm") ) THEN 
     145         DO_2D( 0, 0, 0, 0 ) 
     146            zemx_iwm (ji,jj,1) = 0._wp   ;   zemx_iwm (ji,jj,jpk) = 0._wp 
     147         END_2D 
     148      ENDIF 
     149      IF( iom_use("av_ratio") ) THEN 
     150         DO_2D( 0, 0, 0, 0 ) 
     151            zav_ratio(ji,jj,1) = 0._wp   ;   zav_ratio(ji,jj,jpk) = 0._wp 
     152         END_2D 
     153      ENDIF 
     154      IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN 
     155         DO_2D( 0, 0, 0, 0 ) 
     156            zav_wave (ji,jj,1) = 0._wp   ;   zav_wave (ji,jj,jpk) = 0._wp 
     157         END_2D 
     158      ENDIF 
    143159      ! 
    144160      !                       ! ----------------------------- ! 
    … …  
    148164      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    149165      !                                                 using an exponential decay from the seafloor. 
    150       DO jj = 1, jpj                ! part independent of the level 
    151          DO ji = 1, jpi 
    152             zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    153             zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
    154             IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 
    155          END DO 
    156       END DO 
    157 !!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept_n - sshn 
    158       DO jk = 2, jpkm1              ! complete with the level-dependent part 
    159          DO jj = 1, jpj              
    160             DO ji = 1, jpi 
    161                IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
    162                   zemx_iwm(ji,jj,jk) = 0._wp 
    163                ELSE 
    164                   zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gde3w_n(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )     & 
    165                        &                               - EXP( ( gde3w_n(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) )   & 
    166                        &                            / ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) 
    167                ENDIF 
    168             END DO 
    169          END DO 
    170 !!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point 
     166      DO_2D( 0, 0, 0, 0 ) 
     167         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
     168         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     169         IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj) 
     170      END_2D 
     171!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm) 
     172      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     173         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
     174            zemx_iwm(ji,jj,jk) = 0._wp 
     175         ELSE 
     176            zemx_iwm(ji,jj,jk) = zfact(ji,jj) * (  EXP( ( gde3w(ji,jj,jk  ) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) )     & 
     177                 &                               - EXP( ( gde3w(ji,jj,jk-1) - zhdep(ji,jj) ) / hcri_iwm(ji,jj) ) )   & 
     178                 &                            / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) 
     179         ENDIF 
     180      END_3D 
     181!!gm delta(gde3w) = e3t(:,:,:,Kmm)  !!  Please verify the grid-point position w versus t-point 
    171182!!gm it seems to me that only 1/hcri_iwm  is used ==>  compute it one for all 
    172183 
    173       END DO 
    174184 
    175185      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
    … …  
    180190      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    181191         ! 
    182          zfact(:,:) = 0._wp 
    183          DO jk = 2, jpkm1              ! part independent of the level 
    184             zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    185          END DO 
    186          ! 
    187          DO jj = 1, jpj 
    188             DO ji = 1, jpi 
    189                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    190             END DO 
    191          END DO 
    192          ! 
    193          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    194             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    195          END DO 
     192         DO_2D( 0, 0, 0, 0 ) 
     193            zfact(ji,jj) = 0._wp 
     194         END_2D 
     195         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level 
     196            zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     197         END_3D 
     198         ! 
     199         DO_2D( 0, 0, 0, 0 ) 
     200            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
     201         END_2D 
     202         ! 
     203         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part 
     204            zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     205         END_3D 
    196206         ! 
    197207      CASE ( 2 )               ! Dissipation scales as N^2 
    198208         ! 
    199          zfact(:,:) = 0._wp 
    200          DO jk = 2, jpkm1              ! part independent of the level 
    201             zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    202          END DO 
    203          ! 
    204          DO jj= 1, jpj 
    205             DO ji = 1, jpi 
    206                IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    207             END DO 
    208          END DO 
    209          ! 
    210          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    211             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    212          END DO 
     209         DO_2D( 0, 0, 0, 0 ) 
     210            zfact(ji,jj) = 0._wp 
     211         END_2D 
     212         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level 
     213            zfact(ji,jj) = zfact(ji,jj) + e3w(ji,jj,jk,Kmm) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 
     214         END_3D 
     215         ! 
     216         DO_2D( 0, 0, 0, 0 ) 
     217            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
     218         END_2D 
     219         ! 
     220         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part 
     221            zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zfact(ji,jj) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk) 
     222         END_3D 
    213223         ! 
    214224      END SELECT 
    … …  
    217227      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    218228      ! 
    219       zwkb (:,:,:) = 0._wp 
    220       zfact(:,:)   = 0._wp 
    221       DO jk = 2, jpkm1 
    222          zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    223          zwkb(:,:,jk) = zfact(:,:) 
    224       END DO 
    225 !!gm even better: 
    226 !      DO jk = 2, jpkm1 
    227 !         zwkb(:,:) = zwkb(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) 
    228 !      END DO 
    229 !      zfact(:,:) = zwkb(:,:,jpkm1) 
    230 !!gm or just use zwkb(k=jpk-1) instead of zfact... 
    231 !!gm 
    232       ! 
    233       DO jk = 2, jpkm1 
    234          DO jj = 1, jpj 
    235             DO ji = 1, jpi 
    236                IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
    237                   &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
    238             END DO 
    239          END DO 
    240       END DO 
    241       zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 
    242       ! 
    243       DO jk = 2, jpkm1 
    244          DO jj = 1, jpj 
    245             DO ji = 1, jpi 
    246                IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
    247                   zweight(ji,jj,jk) = 0._wp 
    248                ELSE 
    249                   zweight(ji,jj,jk) = rn2(ji,jj,jk) * hbot_iwm(ji,jj)    & 
    250                      &   * (  EXP( -zwkb(ji,jj,jk) / hbot_iwm(ji,jj) ) - EXP( -zwkb(ji,jj,jk-1) / hbot_iwm(ji,jj) )  ) 
    251                ENDIF 
    252             END DO 
    253          END DO 
    254       END DO 
    255       ! 
    256       zfact(:,:) = 0._wp 
    257       DO jk = 2, jpkm1              ! part independent of the level 
    258          zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    259       END DO 
    260       ! 
    261       DO jj = 1, jpj 
    262          DO ji = 1, jpi 
    263             IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) ) 
    264          END DO 
    265       END DO 
    266       ! 
    267       DO jk = 2, jpkm1              ! complete with the level-dependent part 
    268          zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    269             &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) ) 
    270 !!gm  use of e3t_n just above? 
    271       END DO 
     229      DO_2D( 0, 0, 0, 0 ) 
     230         zwkb(ji,jj,1) = 0._wp 
     231      END_2D 
     232      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     233         zwkb(ji,jj,jk) = zwkb(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(ji,jj,jk) )  ) * wmask(ji,jj,jk) 
     234      END_3D 
     235      DO_2D( 0, 0, 0, 0 ) 
     236         zfact(ji,jj) = zwkb(ji,jj,jpkm1) 
     237      END_2D 
     238      ! 
     239      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     240         IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
     241            &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
     242      END_3D 
     243      DO_2D( 0, 0, 0, 0 ) 
     244         zwkb (ji,jj,1) = zhdep(ji,jj) * wmask(ji,jj,1) 
     245      END_2D 
     246      ! 
     247      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     248         IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization: EXP coast a lot 
     249            zweight(ji,jj,jk) = 0._wp 
     250         ELSE 
     251            zweight(ji,jj,jk) = rn2(ji,jj,jk) * hbot_iwm(ji,jj)    & 
     252               &   * (  EXP( -zwkb(ji,jj,jk) / hbot_iwm(ji,jj) ) - EXP( -zwkb(ji,jj,jk-1) / hbot_iwm(ji,jj) )  ) 
     253         ENDIF 
     254      END_3D 
     255      ! 
     256      DO_2D( 0, 0, 0, 0 ) 
     257         zfact(ji,jj) = 0._wp 
     258      END_2D 
     259      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! part independent of the level 
     260         zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk) 
     261      END_3D 
     262      ! 
     263      DO_2D( 0, 0, 0, 0 ) 
     264         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
     265      END_2D 
     266      ! 
     267      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! complete with the level-dependent part 
     268         zemx_iwm(ji,jj,jk) = zemx_iwm(ji,jj,jk) + zweight(ji,jj,jk) * zfact(ji,jj) * wmask(ji,jj,jk)   & 
     269            &                                                        / ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) 
     270!!gm  use of e3t(ji,jj,:,Kmm) just above? 
     271      END_3D 
    272272      ! 
    273273!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    274274      ! Calculate molecular kinematic viscosity 
    275       znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
    276          &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
    277       DO jk = 2, jpkm1 
    278          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    279       END DO 
     275      DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     276         znu_t(ji,jj,jk) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(ji,jj,jk,jp_tem,Kmm)   & 
     277            &                                     + 0.00694_wp * ts(ji,jj,jk,jp_tem,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)  & 
     278            &                                     + 0.02305_wp * ts(ji,jj,jk,jp_sal,Kmm)  ) * tmask(ji,jj,jk) * r1_rho0 
     279      END_3D 
     280      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     281         znu_w(ji,jj,jk) = 0.5_wp * ( znu_t(ji,jj,jk-1) + znu_t(ji,jj,jk) ) * wmask(ji,jj,jk) 
     282      END_3D 
    280283!!gm end 
    281284      ! 
    282285      ! Calculate turbulence intensity parameter Reb 
    283       DO jk = 2, jpkm1 
    284          zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
    285       END DO 
     286      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     287         zReb(ji,jj,jk) = zemx_iwm(ji,jj,jk) / MAX( 1.e-20_wp, znu_w(ji,jj,jk) * rn2(ji,jj,jk) ) 
     288      END_3D 
    286289      ! 
    287290      ! Define internal wave-induced diffusivity 
    288       DO jk = 2, jpkm1 
    289          zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    290       END DO 
     291      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     292         zav_wave(ji,jj,jk) = znu_w(ji,jj,jk) * zReb(ji,jj,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     293      END_3D 
    291294      ! 
    292295      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    293          DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
    294             DO jj = 1, jpj 
    295                DO ji = 1, jpi 
    296                   IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    297                      zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
    298                   ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 
    299                      zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
    300                   ENDIF 
    301                END DO 
    302             END DO 
    303          END DO 
    304       ENDIF 
    305       ! 
    306       DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    307          zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
    308       END DO 
     296         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     297            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
     298               zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     299            ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 
     300               zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     301            ENDIF 
     302         END_3D 
     303      ENDIF 
     304      ! 
     305      DO_3D( 0, 0, 0, 0, 2, jpkm1 )          ! Bound diffusivity by molecular value and 100 cm2/s 
     306         zav_wave(ji,jj,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(ji,jj,jk) ), 1.e-2_wp  ) * wmask(ji,jj,jk) 
     307      END_3D 
    309308      ! 
    310309      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    311310         zztmp = 0._wp 
    312311!!gm used of glosum 3D.... 
    313          DO jk = 2, jpkm1 
    314             DO jj = 1, jpj 
    315                DO ji = 1, jpi 
    316                   zztmp = zztmp + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   & 
    317                      &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    318                END DO 
    319             END DO 
    320          END DO 
     312         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     313            zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   & 
     314               &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     315         END_3D 
    321316         CALL mpp_sum( 'zdfiwm', zztmp ) 
    322          zztmp = rau0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     317         zztmp = rho0 * zztmp ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
    323318         ! 
    324319         IF(lwp) THEN 
    … …  
    337332      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
    338333         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 
    339          DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb 
    340             DO jj = 1, jpj 
    341                DO ji = 1, jpi 
    342                   ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 
    343                   IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 
    344                      zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) ) 
    345                   ELSE 
    346                      zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk) 
    347                   ENDIF 
    348                END DO 
    349             END DO 
    350          END DO 
     334         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     335            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 
     336            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 
     337               zav_ratio(ji,jj,jk) = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10(ztmp2) - 0.60_wp ) ) 
     338            ELSE 
     339               zav_ratio(ji,jj,jk) = ztmp1 * wmask(ji,jj,jk) 
     340            ENDIF 
     341         END_3D 
    351342         CALL iom_put( "av_ratio", zav_ratio ) 
    352          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
    353             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
    354             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    355             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    356          END DO 
     343         DO_3D( 0, 0, 0, 0, 2, jpkm1 )    !* update momentum & tracer diffusivity with wave-driven mixing 
     344            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) * zav_ratio(ji,jj,jk) 
     345            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     346            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     347         END_3D 
    357348         ! 
    358349      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    359          DO jk = 2, jpkm1 
    360             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) 
    361             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    362             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    363          END DO 
     350         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     351            p_avs(ji,jj,jk) = p_avs(ji,jj,jk) + zav_wave(ji,jj,jk) 
     352            p_avt(ji,jj,jk) = p_avt(ji,jj,jk) + zav_wave(ji,jj,jk) 
     353            p_avm(ji,jj,jk) = p_avm(ji,jj,jk) + zav_wave(ji,jj,jk) 
     354         END_3D 
    364355      ENDIF 
    365356 
    … …  
    368359                                    !* output useful diagnostics: Kz*N^2 ,  
    369360!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 
    370                                     !  vertical integral of rau0 * Kz * N^2 , energy density (zemx_iwm) 
     361                                    !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
    371362      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    372363         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
    373          z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    374          z2d(:,:) = 0._wp 
    375          DO jk = 2, jpkm1 
    376             z2d(:,:) = z2d(:,:) + e3w_n(:,:,jk) * z3d(:,:,jk) * wmask(:,:,jk) 
    377          END DO 
    378          z2d(:,:) = rau0 * z2d(:,:) 
    379          CALL iom_put( "bflx_iwm", z3d ) 
     364         ! Initialisation for iom_put 
     365         DO_2D( 0, 0, 0, 0 ) 
     366            z3d(ji,jj,1) = 0._wp   ;   z3d(ji,jj,jpk) = 0._wp 
     367         END_2D 
     368         z3d(           1:nn_hls,:,:) = 0._wp   ;   z3d(:,           1:nn_hls,:) = 0._wp 
     369         z3d(jpi-nn_hls+1:jpi   ,:,:) = 0._wp   ;   z3d(:,jpj-nn_hls+1:   jpj,:) = 0._wp 
     370         z2d(           1:nn_hls,:  ) = 0._wp   ;   z2d(:,           1:nn_hls  ) = 0._wp 
     371         z2d(jpi-nn_hls+1:jpi   ,:  ) = 0._wp   ;   z2d(:,jpj-nn_hls+1:   jpj  ) = 0._wp 
     372 
     373         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     374            z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) 
     375         END_3D 
     376         DO_2D( 0, 0, 0, 0 ) 
     377            z2d(ji,jj) = 0._wp 
     378         END_2D 
     379         DO_3D( 0, 0, 0, 0, 2, jpkm1 )  
     380            z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk) 
     381         END_3D 
     382         DO_2D( 0, 0, 0, 0 ) 
     383            z2d(ji,jj) = rho0 * z2d(ji,jj) 
     384         END_2D 
     385         CALL iom_put(  "bflx_iwm", z3d ) 
    380386         CALL iom_put( "pcmap_iwm", z2d ) 
    381387         DEALLOCATE( z2d , z3d ) 
    … …  
    383389      CALL iom_put( "emix_iwm", zemx_iwm ) 
    384390       
    385       IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk) 
     391      IF(sn_cfctl%l_prtctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', kdim=jpk) 
    386392      ! 
    387393   END SUBROUTINE zdf_iwm 
    … …  
    414420      !!              de Lavergne et al. in prep., 2017 
    415421      !!---------------------------------------------------------------------- 
    416       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    417       INTEGER  ::   inum         ! local integer 
     422      INTEGER  ::   ifpr               ! dummy loop indices 
     423      INTEGER  ::   inum               ! local integer 
    418424      INTEGER  ::   ios 
    419425      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars 
    420       !! 
    421       NAMELIST/namzdf_iwm/ nn_zpyc, ln_mevar, ln_tsdiff 
    422       !!---------------------------------------------------------------------- 
    423       ! 
    424       REWIND( numnam_ref )              ! Namelist namzdf_iwm in reference namelist : Wave-driven mixing 
     426      ! 
     427      CHARACTER(len=256)            ::   cn_dir                 ! Root directory for location of ssr files 
     428      INTEGER, PARAMETER            ::   jpiwm  = 5             ! maximum number of files to read 
     429      INTEGER, PARAMETER            ::   jp_mpb = 1 
     430      INTEGER, PARAMETER            ::   jp_mpp = 2 
     431      INTEGER, PARAMETER            ::   jp_mpc = 3 
     432      INTEGER, PARAMETER            ::   jp_dsb = 4 
     433      INTEGER, PARAMETER            ::   jp_dsc = 5 
     434      ! 
     435      TYPE(FLD_N), DIMENSION(jpiwm) ::   slf_iwm                ! array of namelist informations 
     436      TYPE(FLD_N)                   ::   sn_mpb, sn_mpp, sn_mpc ! informations about Mixing Power field to be read 
     437      TYPE(FLD_N)                   ::   sn_dsb, sn_dsc         ! informations about Decay Scale field to be read 
     438      TYPE(FLD  ), DIMENSION(jpiwm) ::   sf_iwm                 ! structure of input fields (file informations, fields read) 
     439      ! 
     440      NAMELIST/namzdf_iwm/ nn_zpyc, ln_mevar, ln_tsdiff, & 
     441         &                 cn_dir, sn_mpb, sn_mpp, sn_mpc, sn_dsb, sn_dsc 
     442      !!---------------------------------------------------------------------- 
     443      ! 
    425444      READ  ( numnam_ref, namzdf_iwm, IOSTAT = ios, ERR = 901) 
    426 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist', lwp ) 
    427       ! 
    428       REWIND( numnam_cfg )              ! Namelist namzdf_iwm in configuration namelist : Wave-driven mixing 
     445901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist' ) 
     446      ! 
    429447      READ  ( numnam_cfg, namzdf_iwm, IOSTAT = ios, ERR = 902 ) 
    430 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist', lwp ) 
     448902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist' ) 
    431449      IF(lwm) WRITE ( numond, namzdf_iwm ) 
    432450      ! 
    … …  
    456474      IF( zdf_iwm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' ) 
    457475      ! 
     476      ! store namelist information in an array 
     477      slf_iwm(jp_mpb) = sn_mpb ; slf_iwm(jp_mpp) = sn_mpp ; slf_iwm(jp_mpc) = sn_mpc 
     478      slf_iwm(jp_dsb) = sn_dsb ; slf_iwm(jp_dsc) = sn_dsc 
     479      ! 
     480      DO ifpr= 1, jpiwm 
     481         ALLOCATE( sf_iwm(ifpr)%fnow(jpi,jpj,1)   ) 
     482         IF( slf_iwm(ifpr)%ln_tint )ALLOCATE( sf_iwm(ifpr)%fdta(jpi,jpj,1,2) ) 
     483      END DO 
     484 
     485      ! fill sf_iwm with sf_iwm and control print 
     486      CALL fld_fill( sf_iwm, slf_iwm , cn_dir, 'zdfiwm_init', 'iwm input file', 'namiwm' ) 
     487 
     488      !                             ! hard-coded default definition (to be defined in namelist ?) 
     489      sf_iwm(jp_mpb)%fnow(:,:,1) = 1.e-6 
     490      sf_iwm(jp_mpp)%fnow(:,:,1) = 1.e-6 
     491      sf_iwm(jp_mpc)%fnow(:,:,1) = 1.e-10 
     492      sf_iwm(jp_dsb)%fnow(:,:,1) = 100. 
     493      sf_iwm(jp_dsc)%fnow(:,:,1) = 100. 
     494 
    458495      !                             ! read necessary fields 
    459       CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2] 
    460       CALL iom_get  (inum, jpdom_data, 'field', ebot_iwm, 1 )  
    461       CALL iom_close(inum) 
    462       ! 
    463       CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2] 
    464       CALL iom_get  (inum, jpdom_data, 'field', epyc_iwm, 1 ) 
    465       CALL iom_close(inum) 
    466       ! 
    467       CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2] 
    468       CALL iom_get  (inum, jpdom_data, 'field', ecri_iwm, 1 ) 
    469       CALL iom_close(inum) 
    470       ! 
    471       CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m] 
    472       CALL iom_get  (inum, jpdom_data, 'field', hbot_iwm, 1 ) 
    473       CALL iom_close(inum) 
    474       ! 
    475       CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m] 
    476       CALL iom_get  (inum, jpdom_data, 'field', hcri_iwm, 1 ) 
    477       CALL iom_close(inum) 
    478  
    479       ebot_iwm(:,:) = ebot_iwm(:,:) * ssmask(:,:) 
    480       epyc_iwm(:,:) = epyc_iwm(:,:) * ssmask(:,:) 
    481       ecri_iwm(:,:) = ecri_iwm(:,:) * ssmask(:,:) 
     496      CALL fld_read( nit000, 1, sf_iwm ) 
     497 
     498      ebot_iwm(:,:) = sf_iwm(1)%fnow(:,:,1) * ssmask(:,:) ! energy flux for high-mode wave breaking [W/m2] 
     499      epyc_iwm(:,:) = sf_iwm(2)%fnow(:,:,1) * ssmask(:,:) ! energy flux for pynocline-intensified wave breaking [W/m2] 
     500      ecri_iwm(:,:) = sf_iwm(3)%fnow(:,:,1) * ssmask(:,:) ! energy flux for critical slope wave breaking [W/m2] 
     501      hbot_iwm(:,:) = sf_iwm(4)%fnow(:,:,1)               ! spatially variable decay scale for high-mode wave breaking [m] 
     502      hcri_iwm(:,:) = sf_iwm(5)%fnow(:,:,1)               ! spatially variable decay scale for critical slope wave breaking [m] 
    482503 
    483504      zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) ) 
    484505      zpyc = glob_sum( 'zdfiwm', e1e2t(:,:) * epyc_iwm(:,:) ) 
    485506      zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) ) 
     507 
    486508      IF(lwp) THEN 
    487509         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW' 
Note: See TracChangeset for help on using the changeset viewer.