Changeset 2389


Ignore:
Timestamp:
2010-11-15T16:52:06+01:00 (10 years ago)
Author:
smasson
Message:

ticket #435, bug fix on uslpml and vslpml

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/LDF/ldfslp.F90

    r1515 r2389  
    99   !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90 
    1010   !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates 
     11   !!            3.2  ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML 
    1112   !!---------------------------------------------------------------------- 
    1213#if   defined key_ldfslp   ||   defined key_esopa 
     
    1415   !!   'key_ldfslp'                      Rotation of lateral mixing tensor 
    1516   !!---------------------------------------------------------------------- 
    16    !!   ldf_slp      : compute the slopes of neutral surface 
    17    !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface 
     17   !!   ldf_slp      : compute the slopes of iso-neutral surface 
     18   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface just below the Mixed Layer 
    1819   !!   ldf_slp_init : initialization of the slopes computation 
    1920   !!---------------------------------------------------------------------- 
     
    130131      ENDIF 
    131132       
    132       CALL ldf_slp_mxl( prd, pn2 )                ! Slopes of isopycnal surfaces just below the mixed layer 
    133  
     133 
     134      !                                           !* Local vertical density gradient evaluated from N^2 
     135      DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
     136         zwy(:,:,jk) = zmg * ( prd(:,:,jk) + 1. ) * (    pn2  (:,:,jk) + pn2  (:,:,jk+1)     )   & 
     137            &                                     / MAX( tmask(:,:,jk) + tmask(:,:,jk+1), 1. ) 
     138      END DO 
     139      zwy(:,:,1) = 0.e0                           !  surface value set to zero 
     140 
     141       
     142      CALL ldf_slp_mxl( prd, pn2 )                !* Slopes of isopycnal surfaces just below the mixed layer 
     143       
    134144       
    135145      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd ) 
    136146      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    137147      !                
    138       !                                           !* Local vertical density gradient evaluated from N^2 
    139       DO jk = 2, jpkm1                            !  zwy = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
    140          DO jj = 1, jpj 
    141             DO ji = 1, jpi 
    142                zwy(ji,jj,jk) = zmg * ( prd(ji,jj,jk) + 1. ) * (    pn2  (ji,jj,jk) + pn2  (ji,jj,jk+1)     )   & 
    143                   &                                         / MAX( tmask(ji,jj,jk) + tmask(ji,jj,jk+1), 1. ) 
    144             END DO 
    145          END DO 
    146       END DO 
    147148      DO jk = 2, jpkm1                            !* Slopes at u and v points 
    148149         DO jj = 2, jpjm1 
     
    349350      !!              the mixed layer. 
    350351      !! 
    351       !! ** Method : 
    352       !!      The slope in the i-direction is computed at u- and w-points 
    353       !!   (uslp, wslpi) and the slope in the j-direction is computed at 
    354       !!   v- and w-points (vslp, wslpj). 
    355       !!   They are bounded by 1/100 over the whole ocean, and within the 
    356       !!   surface layer they are bounded by the distance to the surface 
    357       !!   ( slope<= depth/l  where l is the length scale of horizontal 
    358       !!   diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity 
    359       !!   of 10cm/s) 
    360       !! 
    361       !! ** Action :   Compute uslp, wslpi, and vslp, wslpj, the i- and  j-slopes 
    362       !!   of now neutral surfaces at u-, w- and v- w-points, resp. 
    363       !!---------------------------------------------------------------------- 
    364       USE oce , zgru  => ua   ! ua, va used as workspace and set to hor.  
    365       USE oce , zgrv  => va   ! density gradient in ldf_slp 
     352      !! ** Method  :   caution: use zgru, zgrv and zwy computed in ldf_slp 
     353      !! 
     354      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces 
     355      !!                vslpml, wslpjml    just below the mixed layer  
     356      !!                omlmask         :  mixed layer mask 
     357      !!---------------------------------------------------------------------- 
     358      USE oce , zgru  => ua   ! zgru, zgrv (ua, va used as workspace)   
     359      USE oce , zgrv  => va   ! set to i- & j-density gradient in ldf_slp 
     360      USE oce , zwy   => ta   ! set to vertical density gradient at T-point in ldfslp 
    366361      !!  
    367362      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) ::   prd   ! in situ density 
     
    374369      REAL(wp) ::   zau, zbu, zai, zbi   !    -         - 
    375370      REAL(wp) ::   zav, zbv, zaj, zbj   !    -         - 
    376       REAL(wp), DIMENSION(jpi,jpj) ::   zwy   ! 2D workspace 
     371      REAL(wp) ::        zbw             !    -         - 
    377372      !!---------------------------------------------------------------------- 
    378373 
     
    414409      !----------------------------------------------------------------------- 
    415410      ! 
    416       zwy(:,jpj) = 0.e0            !* vertical density gradient for u-slope (from N^2) 
    417       zwy(jpi,:) = 0.e0 
    418 # if defined key_vectopt_loop 
    419       DO jj = 1, 1 
    420          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    421 # else 
    422       DO jj = 1, jpjm1 
    423          DO ji = 1, jpim1 
    424 # endif 
    425             ik = MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) )       ! avoid spurious recirculation  
    426             ik = MIN( ik, jpkm1 )                            ! if ik = jpk take jpkm1 values 
    427             zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   & 
    428                &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. ) 
    429          END DO 
    430       END DO 
    431       CALL lbc_lnk( zwy, 'U', 1. )      ! lateral boundary conditions NO sign change 
    432  
    433411      !                            !* Slope at u points 
    434412# if defined key_vectopt_loop 
     
    443421            ik = MIN( ik, jpkm1 ) 
    444422            zau = 1./ e1u(ji,jj) * zgru(ji,jj,ik) 
    445             zbu = 0.5*( zwy(ji,jj) + zwy(ji+1,jj) ) 
     423            zbu = 0.5*( zwy(ji,jj,ik) + zwy(ji+1,jj,ik) ) 
    446424            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    447425            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     
    452430      END DO 
    453431      CALL lbc_lnk( uslpml, 'U', -1. )      ! lateral boundary conditions (i-gradient => sign change) 
    454  
    455       !                            !* vertical density gradient for v-slope (from N^2) 
    456 # if defined key_vectopt_loop 
    457       DO jj = 1, 1 
    458          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    459 # else 
    460       DO jj = 1, jpjm1 
    461          DO ji = 1, jpim1 
    462 # endif 
    463             ik = MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) 
    464             ik = MIN( ik, jpkm1 ) 
    465             zwy(ji,jj) = zmg * ( prd(ji,jj,ik) + 1. ) * (    pn2  (ji,jj,ik) + pn2  (ji,jj,ik+1)     )   & 
    466                &                                      / MAX( tmask(ji,jj,ik) + tmask(ji,jj,ik+1), 1. ) 
    467          END DO 
    468       END DO 
    469       CALL lbc_lnk( zwy, 'V', 1. )      ! lateral boundary conditions NO sign change 
    470432 
    471433      !                            !* Slope at v points 
     
    481443            ik = MIN( ik,jpkm1 ) 
    482444            zav = 1./ e2v(ji,jj) * zgrv(ji,jj,ik) 
    483             zbv = 0.5*( zwy(ji,jj) + zwy(ji,jj+1) ) 
     445            zbv = 0.5*( zwy(ji,jj,ik) + zwy(ji,jj+1,ik) ) 
    484446            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 
    485447            !                         kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     
    490452      END DO 
    491453      CALL lbc_lnk( vslpml, 'V', -1. )      ! lateral boundary conditions (j-gradient => sign change) 
    492  
    493  
    494       !                            !* vertical density gradient for w-slope (from N^2) 
    495 # if defined key_vectopt_loop 
    496       DO jj = 1, 1 
    497          DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    498 # else 
    499       DO jj = 1, jpj 
    500          DO ji = 1, jpi 
    501 # endif 
    502             ik = nmln(ji,jj) + 1 
    503             ik = MIN( ik, jpk ) 
    504             ikm1 = MAX ( 1, ik-1) 
    505             zwy (ji,jj) = zm05g * pn2 (ji,jj,ik) *     & 
    506                &             ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 
    507          END DO 
    508       END DO 
    509454 
    510455      !                            !* Slopes at w points 
     
    531476            zaj = zcoef2 * (  zgrv(ji,jj  ,ik  ) + zgrv(ji,jj  ,ikm1)   & 
    532477               &            + zgrv(ji,jj-1,ikm1) + zgrv(ji,jj-1,ik  ) ) * tmask (ji,jj,ik) 
     478            ! vertical density gradient at w-point (from N^2) 
     479            zbw = zm05g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 
    533480            ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
    534481            !                   static instability: 
    535482            !                   kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    536             zbi = MIN ( zwy (ji,jj),- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) ) 
    537             zbj = MIN ( zwy (ji,jj), -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) ) 
     483            zbi = MIN ( zbw,- 100.*ABS(zai), -7.e+3/fse3w(ji,jj,ik)*ABS(zai) ) 
     484            zbj = MIN ( zbw, -100.*ABS(zaj), -7.e+3/fse3w(ji,jj,ik)*ABS(zaj) ) 
    538485            ! wslpiml and wslpjml 
    539486            wslpiml (ji,jj) = zai / ( zbi - zeps) * tmask (ji,jj,ik) 
Note: See TracChangeset for help on using the changeset viewer.