Changeset 2401


Ignore:
Timestamp:
2010-11-17T12:09:36+01:00 (10 years ago)
Author:
gm
Message:

v3.3beta: #435: bug correction of the slope just below the ML (uslpiml, vslpjml)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r2371 r2401  
    1010   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates 
    1111   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator 
     12   !!             -   ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML 
    1213   !!---------------------------------------------------------------------- 
    1314#if   defined key_ldfslp   ||   defined key_esopa 
     
    9394      USE oce , zgru  => ua   ! use ua as workspace 
    9495      USE oce , zgrv  => va   ! use va as workspace 
    95       USE oce , zwy   => ta   ! use ta as workspace 
     96      USE oce , zww   => ta   ! use ta as workspace 
    9697      USE oce , zwz   => sa   ! use sa as workspace 
    9798      !! 
     
    103104      INTEGER  ::   ii0, ii1, iku   ! temporary integer 
    104105      INTEGER  ::   ij0, ij1, ikv   ! temporary integer 
    105       REAL(wp) ::   zeps, zmg, zm05g, zalpha        ! temporary scalars 
    106       REAL(wp) ::   zcoef1, zcoef2, zcoef3          !    -         - 
    107       REAL(wp) ::   zcofu , zcofv , zcofw           !    -         - 
    108       REAL(wp) ::   zau, zbu, zai, zbi, z1u, z1wu   !    -         - 
    109       REAL(wp) ::   zav, zbv, zaj, zbj, z1v, z1wv   ! 
    110       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww     ! 3D workspace 
     106      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16     ! local scalars 
     107      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      - 
     108      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
     109      REAL(wp) ::   zck, zfk,      zbw             !   -      - 
     110      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzr   ! 3D workspace 
    111111      !!---------------------------------------------------------------------- 
    112112       
    113       zeps  =  1.e-20                             ! Local constant initialization 
    114       zmg   = -1.0 / grav 
    115       zm05g = -0.5 / grav 
     113      zeps   =  1.e-20_wp        !==   Local constant initialization   ==! 
     114      z1_16  =  1.0_wp / 16._wp 
     115      zm1_g  = -1.0_wp / grav 
     116      zm1_2g = -0.5_wp / grav 
    116117      ! 
    117118      zww(:,:,:) = 0.e0 
    118119      zwz(:,:,:) = 0.e0 
    119       !                                           ! horizontal density gradient computation 
    120       DO jk = 1, jpk 
     120      ! 
     121      DO jk = 1, jpk             !==   i- & j-gradient of density   ==! 
    121122         DO jj = 1, jpjm1 
    122123            DO ji = 1, fs_jpim1   ! vector opt. 
     
    141142         END DO 
    142143      ENDIF 
    143        
    144       CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer 
     144      ! 
     145      zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
     146      DO jk = 2, jpkm1 
     147         !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
     148         !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0 
     149         !                                !    else  tmask(ik+1)  = 0   =>   pn2(ik+1) = 0   =>   zdzr divides by 1 
     150         !                                !          umask(ik+1) /= 0   =>   all pn2  /= 0   =>   zdzr divides by 2 
     151         !                                ! NB: 1/(tmask+1) = (1-.5*tmask)  substitute a / by a *  ==> faster 
     152         zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp )              & 
     153            &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 
     154      END DO 
     155      ! 
     156      !                          !==   Slopes just below the mixed layer   ==! 
     157      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml 
    145158 
    146159       
     
    148161      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    149162      !                
    150       !                                           !* Local vertical density gradient evaluated from N^2 
    151       DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
    152          DO jj = 1, jpj 
    153             DO ji = 1, jpi 
    154                zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * (    pn2  (ji,jj,jk) + pn2  (ji,jj,jk+1)     )   & 
    155                   &                                         / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. ) 
    156             END DO 
    157          END DO 
    158       END DO 
    159163      DO jk = 2, jpkm1                            !* Slopes at u and v points 
    160164         DO jj = 2, jpjm1 
    161165            DO ji = fs_2, fs_jpim1   ! vector opt. 
    162                ! horizontal and vertical density gradient at u- and v-points 
    163                zau = 1. / e1u(ji,jj) * zgru(ji,jj,jk) 
    164                zav = 1. / e2v(ji,jj) * zgrv(ji,jj,jk) 
    165                zbu = 0.5 * ( zwy(ji,jj,jk) + zwy(ji+1,jj  ,jk) ) 
    166                zbv = 0.5 * ( zwy(ji,jj,jk) + zwy(ji  ,jj+1,jk) ) 
    167                ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    168                !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    169                zbu = MIN( zbu, -100.*ABS( zau ), -7.e+3/fse3u(ji,jj,jk)*ABS( zau ) ) 
    170                zbv = MIN( zbv, -100.*ABS( zav ), -7.e+3/fse3v(ji,jj,jk)*ABS( zav ) ) 
    171                ! uslp and vslp output in zwz and zww, resp. 
    172                zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
    173                zwz (ji,jj,jk) = ( ( 1. - zalpha) * zau / ( zbu - zeps )                                           & 
    174                   &                    + zalpha  * uslpml(ji,jj)                                                  & 
    175                   &                              * 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   & 
    176                   &                              / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5. ) ) * umask(ji,jj,jk) 
    177                zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
    178                zww (ji,jj,jk) = ( ( 1. - zalpha) * zav / ( zbv - zeps )                                           & 
    179                   &                    + zalpha  * vslpml(ji,jj)                                                  & 
    180                   &                              * 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   & 
    181                   &                              / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 
     166               !                                      ! horizontal and vertical density gradient at u- and v-points 
     167               zau = zgru(ji,jj,jk) / e1u(ji,jj) 
     168               zav = zgrv(ji,jj,jk) / e2v(ji,jj) 
     169               zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj  ,jk) ) 
     170               zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji  ,jj+1,jk) ) 
     171               !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
     172               !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     173               zbu = MIN(  zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  ) 
     174               zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
     175               !                                      ! uslp and vslp output in zwz and zww, resp. 
     176               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
     177               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
     178               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
     179                  &                   + zfi  * uslpml(ji,jj)                                                     & 
     180                  &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   & 
     181                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 
     182               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
     183                  &                   + zfj  * vslpml(ji,jj)                                                     & 
     184                  &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   & 
     185                  &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 
     186!!gm  modif to suppress omlmask.... (as in Griffies case) 
     187!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0. 
     188!               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 
     189!               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 
     190!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 
     191!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 
     192!               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 
     193!               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 
     194!!gm end modif 
    182195            END DO 
    183196         END DO 
     
    185198      CALL lbc_lnk( zwz, 'U', -1. )   ;   CALL lbc_lnk( zww, 'V', -1. )      ! lateral boundary conditions 
    186199      ! 
    187       zcofu = 1. / 16.                            !* horizontal Shapiro filter 
    188       zcofv = 1. / 16. 
     200      !                                            !* horizontal Shapiro filter 
    189201      DO jk = 2, jpkm1 
    190202         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 only 
    191203            DO ji = 2, jpim1   
    192                uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
     204               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
    193205                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
    194206                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      & 
    195207                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    & 
    196208                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       ) 
    197                vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
     209               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
    198210                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      & 
    199211                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      & 
     
    204216         DO jj = 3, jpj-2                               ! other rows 
    205217            DO ji = fs_2, fs_jpim1   ! vector opt. 
    206                uslp(ji,jj,jk) = zcofu * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
     218               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
    207219                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
    208220                  &                       + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)      & 
    209221                  &                       +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )    & 
    210222                  &                       + 4.*  zwz(ji  ,jj  ,jk)                       ) 
    211                vslp(ji,jj,jk) = zcofv * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
     223               vslp(ji,jj,jk) = z1_16 * (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)      & 
    212224                  &                       +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)      & 
    213225                  &                       + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)      & 
     
    219231         DO jj = 2, jpjm1 
    220232            DO ji = fs_2, fs_jpim1   ! vector opt. 
    221                z1u  = ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk) )*.5 
    222                z1v  = ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk) )*.5 
    223                z1wu = ( umask(ji,jj,jk)   + umask(ji,jj,jk+1) )*.5 
    224                z1wv = ( vmask(ji,jj,jk)   + vmask(ji,jj,jk+1) )*.5 
    225                uslp(ji,jj,jk) = uslp(ji,jj,jk) * z1u * z1wu 
    226                vslp(ji,jj,jk) = vslp(ji,jj,jk) * z1v * z1wv 
     233               uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
     234                  &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp 
     235               vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
     236                  &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp 
    227237            END DO 
    228238         END DO 
     
    233243      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd ) 
    234244      !                
    235       !                                           !* Local vertical density gradient evaluated from N^2 
    236       DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - mk ( prd ) / grav * pn2 -- at w point 
    237          DO jj = 1, jpj 
    238             DO ji = 1, jpi 
    239                zwy(ji,jj,jk) = zm05g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 
    240             END DO 
    241          END DO 
    242       END DO 
    243       DO jk = 2, jpkm1                            !* Slopes at w point 
     245      DO jk = 2, jpkm1 
    244246         DO jj = 2, jpjm1 
    245247            DO ji = fs_2, fs_jpim1   ! vector opt. 
    246                !                                        ! horizontal density i-gradient at w-points 
    247                zcoef1 = MAX( zeps, umask(ji-1,jj,jk  )+umask(ji,jj,jk  )    & 
    248                   &               +umask(ji-1,jj,jk-1)+umask(ji,jj,jk-1) ) 
    249                zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) ) 
    250                zai = zcoef1 * (  zgru(ji  ,jj,jk  ) + zgru(ji  ,jj,jk-1)   & 
    251                   &            + zgru(ji-1,jj,jk-1) + zgru(ji-1,jj,jk  ) ) * tmask (ji,jj,jk) 
    252                !                                        ! horizontal density j-gradient at w-points 
    253                zcoef2 = MAX( zeps, vmask(ji,jj-1,jk  )+vmask(ji,jj,jk-1)   & 
    254                   &               +vmask(ji,jj-1,jk-1)+vmask(ji,jj,jk  ) ) 
    255                zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) ) 
    256                zaj = zcoef2 * (  zgrv(ji,jj  ,jk  ) + zgrv(ji,jj  ,jk-1)   & 
    257                   &            + zgrv(ji,jj-1,jk-1) + zgrv(ji,jj-1,jk  ) ) * tmask (ji,jj,jk) 
     248               !                                  !* Local vertical density gradient evaluated from N^2 
     249               zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 
     250               !                                  !* Slopes at w point 
     251               !                                        ! i- & j-gradient of density at w-points 
     252               zci = MAX(  umask(ji-1,jj,jk  ) + umask(ji,jj,jk  )           & 
     253                  &      + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps  ) * e1t(ji,jj) 
     254               zcj = MAX(  vmask(ji,jj-1,jk  ) + vmask(ji,jj,jk-1)           & 
     255                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj) 
     256               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           & 
     257                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk) 
     258               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           & 
     259                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk) 
    258260               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
    259                !                                        ! static instability: kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    260                zbi = MIN( zwy (ji,jj,jk),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,jk)*ABS(zai) ) 
    261                zbj = MIN( zwy (ji,jj,jk), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,jk)*ABS(zaj) ) 
    262                !                                        ! wslpi and wslpj output in zwz and zww, resp. 
    263                zalpha = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) 
    264                zcoef3 = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 
    265                zwz(ji,jj,jk) = (     zai / ( zbi - zeps)  * ( 1. - zalpha )   & 
    266                   &             + zcoef3 * wslpiml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk) 
    267                zww(ji,jj,jk) = (     zaj / ( zbj - zeps)  * ( 1. - zalpha )   & 
    268                   &             + zcoef3 * wslpjml(ji,jj) *        zalpha   ) * tmask (ji,jj,jk) 
    269             END DO 
    270          END DO 
    271       END DO 
    272       CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions on zwz and zww 
     261               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     262               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  ) 
     263               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  ) 
     264               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
     265               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
     266               zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 
     267               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
     268               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
     269 
     270!!gm  modif to suppress omlmask....  (as in Griffies operator) 
     271!               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0. 
     272!               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 
     273!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. ) 
     274!               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 
     275!               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 
     276!!gm end modif 
     277            END DO 
     278         END DO 
     279      END DO 
     280      CALL lbc_lnk( zwz, 'T', -1. )   ;    CALL lbc_lnk( zww, 'T', -1. )      ! lateral boundary conditions 
    273281      ! 
    274282      !                                           !* horizontal Shapiro filter 
     
    276284         DO jj = 2, jpjm1, jpj-3                        ! rows jj=2 and =jpjm1 
    277285            DO ji = 2, jpim1 
    278                zcofw = tmask(ji,jj,jk) / 16. 
    279286               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
    280287                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
    281288                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
    282289                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
    283                   &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw 
     290                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
    284291 
    285292               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
     
    287294                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
    288295                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
    289                   &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw 
     296                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
    290297            END DO 
    291298         END DO   
    292299         DO jj = 3, jpj-2                               ! other rows 
    293300            DO ji = fs_2, fs_jpim1   ! vector opt. 
    294                zcofw = tmask(ji,jj,jk) / 16. 
    295301               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
    296302                  &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
    297303                  &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
    298304                  &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
    299                   &                + 4.*  zwz(ji  ,jj  ,jk)                        ) * zcofw 
     305                  &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
    300306 
    301307               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
     
    303309                  &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
    304310                  &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
    305                   &                + 4.*  zww(ji  ,jj  ,jk)                        ) * zcofw 
     311                  &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
    306312            END DO 
    307313         END DO 
     
    309315         DO jj = 2, jpjm1 
    310316            DO ji = fs_2, fs_jpim1   ! vector opt. 
    311                z1u = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) *.5 
    312                z1v = ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) *.5 
    313                wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * z1u * z1v 
    314                wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * z1u * z1v 
     317               zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
     318                  &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
     319               wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 
     320               wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 
    315321            END DO 
    316322         END DO 
     
    573579 
    574580 
    575    SUBROUTINE ldf_slp_mxl( prd, pn2 ) 
     581   SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr ) 
    576582      !!---------------------------------------------------------------------- 
    577583      !!                  ***  ROUTINE ldf_slp_mxl  *** 
     
    580586      !!              the mixed layer. 
    581587      !! 
    582       !! ** Method : 
    583       !!      The slope in the i-direction is computed at u- and w-points 
    584       !!   (uslp, wslpi) and the slope in the j-direction is computed at 
    585       !!   v- and w-points (vslp, wslpj). 
    586       !!   They are bounded by 1/100 over the whole ocean, and within the 
    587       !!   surface layer they are bounded by the distance to the surface 
    588       !!   ( slope<= depth/l  where l is the length scale of horizontal 
    589       !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity 
    590       !!   of 10cm/s) 
    591       !! 
    592       !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes 
    593       !!   of now neutral surfaces at u-, w- and v- w-points, resp. 
    594       !!---------------------------------------------------------------------- 
    595       USE oce , zgru  => ua   ! ua, va used as workspace and set to hor.  
    596       USE oce , zgrv  => va   ! density gradient in ldf_slp 
    597       !!  
    598       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density 
    599       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2   ! Brunt-Vaisala frequency (locally ref.) 
    600       !! 
    601       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    602       INTEGER  ::   ik, ikm1             ! temporary integers 
    603       REAL(wp) ::   zeps, zmg, zm05g     ! temporary scalars 
    604       REAL(wp) ::   zcoef1, zcoef2       !    -         - 
    605       REAL(wp) ::   zau, zbu, zai, zbi   !    -         - 
    606       REAL(wp) ::   zav, zbv, zaj, zbj   !    -         - 
    607       REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace 
    608       !!---------------------------------------------------------------------- 
    609  
    610       zeps  =  1.e-20                    ! Local constant initialization 
    611       zmg   = -1.0 / grav 
    612       zm05g = -0.5 / grav 
     588      !! ** Method  :   The slope in the i-direction is computed at u- & w-points 
     589      !!              (uslpml, wslpiml) and the slope in the j-direction is computed 
     590      !!              at v- and w-points (vslpml, wslpjml) with the same bounds as 
     591      !!              in ldf_slp. 
     592      !! 
     593      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces 
     594      !!                vslpml, wslpjml    just below the mixed layer  
     595      !!                omlmask         :  mixed layer mask 
     596      !!---------------------------------------------------------------------- 
     597      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd            ! in situ density 
     598      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   pn2            ! Brunt-Vaisala frequency (locally ref.) 
     599      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   p_gru, p_grv   ! i- & j-gradient of density (u- & v-pts) 
     600      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   p_dzr          ! z-gradient of density      (T-point) 
     601      !! 
     602      INTEGER  ::   ji , jj , jk         ! dummy loop indices 
     603      INTEGER  ::   iku, ikv, ik, ikm1   ! local integers 
     604      REAL(wp) ::   zeps, zm1_g, zm1_2g            ! local scalars 
     605      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      - 
     606      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
     607      REAL(wp) ::   zck, zfk,      zbw             !   -      - 
     608      !!---------------------------------------------------------------------- 
     609 
     610      zeps   =  1.e-20_wp        !==   Local constant initialization   ==! 
     611      zm1_g  = -1.0_wp / grav 
     612      zm1_2g = -0.5_wp / grav 
    613613      ! 
    614614      uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0 
     
    616616      wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0 
    617617      wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0 
    618  
    619       !                                  ! surface mixed layer mask 
    620       DO jk = 1, jpk                          ! =1 inside the mixed layer, =0 otherwise 
     618      ! 
     619      !                          !==   surface mixed layer mask   ! 
     620      DO jk = 1, jpk                      ! =1 inside the mixed layer, =0 otherwise 
    621621# if defined key_vectopt_loop 
    622622         DO jj = 1, 1 
     
    645645      !----------------------------------------------------------------------- 
    646646      ! 
    647       zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2) 
    648       zwy(jpi,:) = 0.e0 
    649 # if defined key_vectopt_loop 
    650       DO jj = 1, 1 
    651          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    652 # else 
    653       DO jj = 1, jpjm1 
    654          DO ji = 1, jpim1 
    655 # endif 
    656             ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )       ! avoid spurious recirculation  
    657             ik = MIN( ik, jpkm1 )                            ! if ik = jpk take jpkm1 values 
    658             zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   & 
    659                &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. ) 
    660          END DO 
    661       END DO 
    662       CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions NO sign change 
    663  
    664       !                            !* Slope at u points 
    665647# if defined key_vectopt_loop 
    666648      DO jj = 1, 1 
     
    670652         DO ji = 2, jpim1 
    671653# endif 
    672             ! horizontal and vertical density gradient at u-points 
    673             ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) 
    674             ik = MIN( ik, jpkm1 ) 
    675             zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik) 
    676             zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) ) 
    677             ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    678             !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    679             zbu = MIN( zbu, -100.*ABS(zau), -7.e+3/fse3u(ji,jj,ik)*ABS(zau) ) 
    680             ! uslpml 
    681             uslpml (ji,jj) = zau / ( zbu - zeps ) * umask (ji,jj,ik) 
    682          END DO 
    683       END DO 
    684       CALL lbc_lnk( uslpml, 'U', -1. )      ! lateral boundary conditions (i-gradient => sign change) 
    685  
    686       !                            !* vertical density gradient for v-slope (from N^2) 
    687 # if defined key_vectopt_loop 
    688       DO jj = 1, 1 
    689          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    690 # else 
    691       DO jj = 1, jpjm1 
    692          DO ji = 1, jpim1 
    693 # endif 
    694             ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) 
    695             ik = MIN( ik, jpkm1 ) 
    696             zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   & 
    697                &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. ) 
    698          END DO 
    699       END DO 
    700       CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions NO sign change 
    701  
    702       !                            !* Slope at v points 
    703 # if defined key_vectopt_loop 
    704       DO jj = 1, 1 
    705          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    706 # else 
    707       DO jj = 2, jpjm1 
    708          DO ji = 2, jpim1 
    709 # endif 
    710             ! horizontal and vertical density gradient at v-points 
    711             ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) 
    712             ik = MIN( ik,jpkm1 ) 
    713             zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik) 
    714             zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) ) 
    715             ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    716             !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    717             zbv = MIN( zbv, -100.*ABS(zav), -7.e+3/fse3v(ji,jj,ik)*ABS( zav ) ) 
    718             ! vslpml 
    719             vslpml (ji,jj) = zav / ( zbv - zeps ) * vmask (ji,jj,ik) 
    720          END DO 
    721       END DO 
    722       CALL lbc_lnk( vslpml, 'V', -1. )      ! lateral boundary conditions (j-gradient => sign change) 
    723  
    724  
    725       !                            !* vertical density gradient for w-slope (from N^2) 
    726 # if defined key_vectopt_loop 
    727       DO jj = 1, 1 
    728          DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    729 # else 
    730       DO jj = 1, jpj 
    731          DO ji = 1, jpi 
    732 # endif 
    733             ik = nmln(ji,jj) + 1 
    734             ik = MIN( ik, jpk ) 
    735             ikm1 = MAX ( 1, ik-1) 
    736             zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     & 
    737                &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 
    738          END DO 
    739       END DO 
    740  
    741       !                            !* Slopes at w points 
    742 # if defined key_vectopt_loop 
    743       DO jj = 1, 1 
    744          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    745 # else 
    746       DO jj = 2, jpjm1 
    747          DO ji = 2, jpim1 
    748 # endif 
    749             ik = nmln(ji,jj) + 1 
    750             ik = MIN( ik, jpk ) 
    751             ikm1 = MAX ( 1, ik-1 ) 
    752             ! horizontal density i-gradient at w-points 
    753             zcoef1 = MAX( zeps, umask(ji-1,jj,ik  )+umask(ji,jj,ik  )   & 
    754                &               +umask(ji-1,jj,ikm1)+umask(ji,jj,ikm1) ) 
    755             zcoef1 = 1. / ( zcoef1 * e1t (ji,jj) ) 
    756             zai = zcoef1 * (  zgru(ji  ,jj,ik  ) + zgru(ji  ,jj,ikm1)   & 
    757                &            + zgru(ji-1,jj,ikm1) + zgru(ji-1,jj,ik  ) ) * tmask (ji,jj,ik) 
    758             ! horizontal density j-gradient at w-points 
    759             zcoef2 = MAX( zeps, vmask(ji,jj-1,ik  )+vmask(ji,jj,ikm1)    & 
    760                &               +vmask(ji,jj-1,ikm1)+vmask(ji,jj,ik  ) ) 
    761             zcoef2 = 1.0 / ( zcoef2 *  e2t (ji,jj) ) 
    762             zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   & 
    763                &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik) 
    764             ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
    765             !                   static instability: 
    766             !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    767             zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) ) 
    768             zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) ) 
    769             ! wslpiml and wslpjml 
    770             wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik) 
    771             wslpjml (ji,jj) = zaj / ( zbj - zeps) * tmask (ji,jj,ik) 
    772          END DO 
    773       END DO 
    774       CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )      ! lateral boundary conditions 
     654            !                    !==   Slope at u- & v-points just below the Mixed Layer   ==! 
     655            ! 
     656            !                          !- vertical density gradient for u- and v-slopes (from dzr at T-point) 
     657            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
     658            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !  
     659            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) ) 
     660            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) ) 
     661            !                          !- horizontal density gradient at u- & v-points 
     662            zau = p_gru(ji,jj,iku) / e1u(ji,jj) 
     663            zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 
     664            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 
     665            !                                kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     666            zbu = MIN(  zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,ik)* ABS( zau )  ) 
     667            zbv = MIN(  zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ik)* ABS( zav )  ) 
     668            !                          !- Slope at u- & v-points (uslpml, vslpml) 
     669            uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,ik) 
     670            vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ik) 
     671            ! 
     672            !                    !==   i- & j-slopes at w-points just below the Mixed Layer   ==! 
     673            ! 
     674            ik   = MIN( nmln(ji,jj) + 1, jpk ) 
     675            ikm1 = MAX( 1, ik-1 ) 
     676            !                          !- vertical density gradient for w-slope (from N^2) 
     677            zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 
     678            !                          !- horizontal density i- & j-gradient at w-points 
     679            zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           & 
     680               &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj)  
     681            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           & 
     682               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj) 
     683            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)           & 
     684               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask(ji,jj,ik) 
     685            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           & 
     686               &       + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1)  ) / zcj  * tmask(ji,jj,ik) 
     687            !                          !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
     688            !                             kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     689            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai )  ) 
     690            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  ) 
     691            !                          !- i- & j-slope at w-points (wslpiml, wslpjml) 
     692            wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 
     693            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 
     694         END DO 
     695      END DO 
     696!!gm this lbc_lnk should be useless.... 
     697      CALL lbc_lnk( uslpml , 'U', -1. )   ;   CALL lbc_lnk( vslpml , 'V', -1. )   ! lateral boundary cond. (sign change) 
     698      CALL lbc_lnk( wslpiml, 'W', -1. )   ;   CALL lbc_lnk( wslpjml, 'W', -1. )   ! lateral boundary conditions 
    775699      ! 
    776700   END SUBROUTINE ldf_slp_mxl 
Note: See TracChangeset for help on using the changeset viewer.