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 13540 for NEMO/branches/2020/r12377_ticket2386/src/OCE/ZDF/zdfiwm.F90 – NEMO

Ignore:
Timestamp:
2020-09-29T12:41:06+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2386: update to latest trunk

Location:
NEMO/branches/2020/r12377_ticket2386
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r12377_ticket2386

    • 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 
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@13507        sette 
  • NEMO/branches/2020/r12377_ticket2386/src/OCE/ZDF/zdfiwm.F90

    r12511 r13540  
    2323   USE phycst         ! physical constants 
    2424   ! 
     25   USE fldread        ! field read 
    2526   USE prtctl         ! Print control 
    2627   USE in_out_manager ! I/O manager 
     
    5051   !! * Substitutions 
    5152#  include "do_loop_substitute.h90" 
     53#  include "domzgr_substitute.h90" 
    5254   !!---------------------------------------------------------------------- 
    5355   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    9395      !!                 2. Pycnocline-intensified low-mode dissipation 
    9496      !!                     zemx_iwm(z) = ( epyc_iwm / rho0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
    95       !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
     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 
     
    99101      !!                 3. WKB-height dependent high mode dissipation 
    100102      !!                     zemx_iwm(z) = ( ebot_iwm / rho0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm) 
    101       !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) ) 
     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:  
     
    138140      !!---------------------------------------------------------------------- 
    139141      ! 
    140       !                       !* Set to zero the 1st and last vertical levels of appropriate variables 
    141       zemx_iwm (:,:,1) = 0._wp   ;   zemx_iwm (:,:,jpk) = 0._wp 
    142       zav_ratio(:,:,1) = 0._wp   ;   zav_ratio(:,:,jpk) = 0._wp 
    143       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 
    144159      ! 
    145160      !                       ! ----------------------------- ! 
     
    149164      !                       !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
    150165      !                                                 using an exponential decay from the seafloor. 
    151       DO_2D_11_11 
     166      DO_2D( 0, 0, 0, 0 )             ! part independent of the level 
    152167         zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    153168         zfact(ji,jj) = rho0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  ) 
     
    155170      END_2D 
    156171!!gm gde3w ==>>>  check for ssh taken into account.... seem OK gde3w_n=gdept(:,:,:,Kmm) - ssh(:,:,Kmm) 
    157       DO_3D_11_11( 2, jpkm1 ) 
     172      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! complete with the level-dependent part 
    158173         IF ( zfact(ji,jj) == 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
    159174            zemx_iwm(ji,jj,jk) = 0._wp 
     
    175190      CASE ( 1 )               ! Dissipation scales as N (recommended) 
    176191         ! 
    177          zfact(:,:) = 0._wp 
    178          DO jk = 2, jpkm1              ! part independent of the level 
    179             zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    180          END DO 
    181          ! 
    182          DO_2D_11_11 
     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 ) 
    183200            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    184201         END_2D 
    185202         ! 
    186          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    187             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    188          END DO 
     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 
    189206         ! 
    190207      CASE ( 2 )               ! Dissipation scales as N^2 
    191208         ! 
    192          zfact(:,:) = 0._wp 
    193          DO jk = 2, jpkm1              ! part independent of the level 
    194             zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    195          END DO 
    196          ! 
    197          DO_2D_11_11 
     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 ) 
    198217            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    199218         END_2D 
    200219         ! 
    201          DO jk = 2, jpkm1              ! complete with the level-dependent part 
    202             zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
    203          END DO 
     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 
    204223         ! 
    205224      END SELECT 
     
    208227      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
    209228      ! 
    210       zwkb (:,:,:) = 0._wp 
    211       zfact(:,:)   = 0._wp 
    212       DO jk = 2, jpkm1 
    213          zfact(:,:) = zfact(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
    214          zwkb(:,:,jk) = zfact(:,:) 
    215       END DO 
    216 !!gm even better: 
    217 !      DO jk = 2, jpkm1 
    218 !         zwkb(:,:) = zwkb(:,:) + e3w(:,:,jk,Kmm) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) 
    219 !      END DO 
    220 !      zfact(:,:) = zwkb(:,:,jpkm1) 
    221 !!gm or just use zwkb(k=jpk-1) instead of zfact... 
    222 !!gm 
    223       ! 
    224       DO_3D_11_11( 2, jpkm1 ) 
     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 ) 
    225240         IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
    226241            &                                     * wmask(ji,jj,jk) / zfact(ji,jj) 
    227242      END_3D 
    228       zwkb(:,:,1) = zhdep(:,:) * wmask(:,:,1) 
    229       ! 
    230       DO_3D_11_11( 2, jpkm1 ) 
    231          IF ( rn2(ji,jj,jk) <= 0._wp .OR. wmask(ji,jj,jk) == 0._wp ) THEN   ! optimization 
     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 
    232249            zweight(ji,jj,jk) = 0._wp 
    233250         ELSE 
     
    237254      END_3D 
    238255      ! 
    239       zfact(:,:) = 0._wp 
    240       DO jk = 2, jpkm1              ! part independent of the level 
    241          zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
    242       END DO 
    243       ! 
    244       DO_2D_11_11 
     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 ) 
    245264         IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rho0 * zfact(ji,jj) ) 
    246265      END_2D 
    247266      ! 
    248       DO jk = 2, jpkm1              ! complete with the level-dependent part 
    249          zemx_iwm(:,:,jk) = zemx_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
    250             &                                / ( gde3w(:,:,jk) - gde3w(:,:,jk-1) ) 
    251 !!gm  use of e3t(:,:,:,Kmm) just above? 
    252       END DO 
     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 
    253272      ! 
    254273!!gm  this is to be replaced by just a constant value znu=1.e-6 m2/s 
    255274      ! Calculate molecular kinematic viscosity 
    256       znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * ts(:,:,:,jp_tem,Kmm) + 0.00694_wp * ts(:,:,:,jp_tem,Kmm) * ts(:,:,:,jp_tem,Kmm)  & 
    257          &                                  + 0.02305_wp * ts(:,:,:,jp_sal,Kmm)  ) * tmask(:,:,:) * r1_rho0 
    258       DO jk = 2, jpkm1 
    259          znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
    260       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 
    261283!!gm end 
    262284      ! 
    263285      ! Calculate turbulence intensity parameter Reb 
    264       DO jk = 2, jpkm1 
    265          zReb(:,:,jk) = zemx_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
    266       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 
    267289      ! 
    268290      ! Define internal wave-induced diffusivity 
    269       DO jk = 2, jpkm1 
    270          zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
    271       END DO 
    272       ! 
    273       IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
    274          DO_3D_11_11( 2, jpkm1 ) 
     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 
     294      ! 
     295      IF( ln_mevar ) THEN                ! Variable mixing efficiency case : modify zav_wave in the 
     296         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
    275297            IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
    276298               zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     
    281303      ENDIF 
    282304      ! 
    283       DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
    284          zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
    285       END DO 
     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 
    286308      ! 
    287309      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
    288310         zztmp = 0._wp 
    289311!!gm used of glosum 3D.... 
    290          DO_3D_11_11( 2, jpkm1 ) 
     312         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    291313            zztmp = zztmp + e3w(ji,jj,jk,Kmm) * e1e2t(ji,jj)   & 
    292314               &          * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     
    308330      !                          ! ----------------------- ! 
    309331      !       
    310       IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
     332      IF( ln_tsdiff ) THEN                !* Option for differential mixing of salinity and temperature 
    311333         ztmp1 = 0.505_wp + 0.495_wp * TANH( 0.92_wp * ( LOG10( 1.e-20_wp ) - 0.60_wp ) ) 
    312          DO_3D_11_11( 2, jpkm1 ) 
     334         DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! Calculate S/T diffusivity ratio as a function of Reb 
    313335            ztmp2 = zReb(ji,jj,jk) * 5._wp * r1_6 
    314336            IF ( ztmp2 > 1.e-20_wp .AND. wmask(ji,jj,jk) == 1._wp ) THEN 
     
    319341         END_3D 
    320342         CALL iom_put( "av_ratio", zav_ratio ) 
    321          DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
    322             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
    323             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    324             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    325          END DO 
    326          ! 
    327       ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
    328          DO jk = 2, jpkm1 
    329             p_avs(:,:,jk) = p_avs(:,:,jk) + zav_wave(:,:,jk) 
    330             p_avt(:,:,jk) = p_avt(:,:,jk) + zav_wave(:,:,jk) 
    331             p_avm(:,:,jk) = p_avm(:,:,jk) + zav_wave(:,:,jk) 
    332          END DO 
    333       ENDIF 
    334  
    335       !                             !* output internal wave-driven mixing coefficient 
     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 
     348         ! 
     349      ELSE                                !* update momentum & tracer diffusivity with wave-driven mixing 
     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 
     355      ENDIF 
     356 
     357      !                                   !* output internal wave-driven mixing coefficient 
    336358      CALL iom_put( "av_wave", zav_wave ) 
    337                                     !* output useful diagnostics: Kz*N^2 ,  
     359                                          !* output useful diagnostics: Kz*N^2 ,  
    338360!!gm Kz*N2 should take into account the ratio avs/avt if it is used.... (see diaar5) 
    339                                     !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
     361                                          !  vertical integral of rho0 * Kz * N^2 , energy density (zemx_iwm) 
    340362      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN 
    341363         ALLOCATE( z2d(jpi,jpj) , z3d(jpi,jpj,jpk) ) 
    342          z3d(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
    343          z2d(:,:) = 0._wp 
    344          DO jk = 2, jpkm1 
    345             z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk) 
    346          END DO 
    347          z2d(:,:) = rho0 * z2d(:,:) 
    348          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 ) 
    349386         CALL iom_put( "pcmap_iwm", z2d ) 
    350387         DEALLOCATE( z2d , z3d ) 
     
    383420      !!              de Lavergne et al. in prep., 2017 
    384421      !!---------------------------------------------------------------------- 
    385       INTEGER  ::   inum         ! local integer 
     422      INTEGER  ::   ifpr               ! dummy loop indices 
     423      INTEGER  ::   inum               ! local integer 
    386424      INTEGER  ::   ios 
    387425      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars 
    388       !! 
    389       NAMELIST/namzdf_iwm/ nn_zpyc, ln_mevar, ln_tsdiff 
     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 
    390442      !!---------------------------------------------------------------------- 
    391443      ! 
     
    422474      IF( zdf_iwm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' ) 
    423475      ! 
     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 
    424495      !                             ! read necessary fields 
    425       CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2] 
    426       CALL iom_get  (inum, jpdom_data, 'field', ebot_iwm, 1 )  
    427       CALL iom_close(inum) 
    428       ! 
    429       CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2] 
    430       CALL iom_get  (inum, jpdom_data, 'field', epyc_iwm, 1 ) 
    431       CALL iom_close(inum) 
    432       ! 
    433       CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2] 
    434       CALL iom_get  (inum, jpdom_data, 'field', ecri_iwm, 1 ) 
    435       CALL iom_close(inum) 
    436       ! 
    437       CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m] 
    438       CALL iom_get  (inum, jpdom_data, 'field', hbot_iwm, 1 ) 
    439       CALL iom_close(inum) 
    440       ! 
    441       CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m] 
    442       CALL iom_get  (inum, jpdom_data, 'field', hcri_iwm, 1 ) 
    443       CALL iom_close(inum) 
    444  
    445       ebot_iwm(:,:) = ebot_iwm(:,:) * ssmask(:,:) 
    446       epyc_iwm(:,:) = epyc_iwm(:,:) * ssmask(:,:) 
    447       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] 
    448503 
    449504      zbot = glob_sum( 'zdfiwm', e1e2t(:,:) * ebot_iwm(:,:) ) 
    450505      zpyc = glob_sum( 'zdfiwm', e1e2t(:,:) * epyc_iwm(:,:) ) 
    451506      zcri = glob_sum( 'zdfiwm', e1e2t(:,:) * ecri_iwm(:,:) ) 
     507 
    452508      IF(lwp) THEN 
    453509         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW' 
Note: See TracChangeset for help on using the changeset viewer.