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 5581 for branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 – NEMO

Ignore:
Timestamp:
2015-07-10T13:28:53+02:00 (9 years ago)
Author:
timgraham
Message:

Merged head of trunk into branch

File:
1 edited

Legend:

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

    r4488 r5581  
    2828   USE zdfmxl         ! mixed layer depth 
    2929   USE eosbn2         ! equation of states 
     30   ! 
     31   USE in_out_manager ! I/O manager 
    3032   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    31    USE in_out_manager ! I/O manager 
    3233   USE prtctl         ! Print control 
    3334   USE wrk_nemo       ! work arrays 
     
    108109      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
    109110      REAL(wp) ::   zck, zfk,      zbw             !   -      - 
     111      REAL(wp) ::   zdepv, zdepu         !   -      - 
    110112      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww 
    111113      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr 
    112114      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zgru, zgrv 
     115      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhmlpu, zhmlpv 
    113116      !!---------------------------------------------------------------------- 
    114117      ! 
     
    116119      ! 
    117120      CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 
     121      CALL wrk_alloc( jpi,jpj, zhmlpu, zhmlpv ) 
    118122 
    119123      IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN  
     
    136140         END DO 
    137141         IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level 
    138 # if defined key_vectopt_loop 
    139             DO jj = 1, 1 
    140                DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    141 # else 
    142142            DO jj = 1, jpjm1 
    143143               DO ji = 1, jpim1 
    144 # endif 
    145144                  zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
    146145                  zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     
    148147            END DO 
    149148         ENDIF 
    150          ! 
    151          zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
     149         IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level 
     150            DO jj = 1, jpjm1 
     151               DO ji = 1, jpim1 
     152               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
     153               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 
     154               END DO 
     155            END DO 
     156         ENDIF 
     157         ! 
     158         !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
     159         ! interior value 
    152160         DO jk = 2, jpkm1 
    153161            !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
     
    159167               &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 
    160168         END DO 
     169         ! surface initialisation  
     170         zdzr(:,:,1) = 0._wp  
     171         IF ( ln_isfcav ) THEN 
     172            ! if isf need to overwrite the interior value at at the first ocean point 
     173            DO jj = 1, jpjm1 
     174               DO ji = 1, jpim1 
     175                  zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 
     176               END DO 
     177            END DO 
     178         END IF 
    161179         ! 
    162180         !                          !==   Slopes just below the mixed layer   ==! 
     
    167185         ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    168186         ! 
     187         IF ( ln_isfcav ) THEN 
     188            DO jj = 2, jpjm1 
     189               DO ji = fs_2, fs_jpim1   ! vector opt. 
     190                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
     191                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp) 
     192                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp) 
     193                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
     194                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp) 
     195                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp) 
     196               ENDDO 
     197            ENDDO 
     198         ELSE 
     199            DO jj = 2, jpjm1 
     200               DO ji = fs_2, fs_jpim1   ! vector opt. 
     201                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) 
     202                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) 
     203               ENDDO 
     204            ENDDO 
     205         END IF 
    169206         DO jk = 2, jpkm1                            !* Slopes at u and v points 
    170207            DO jj = 2, jpjm1 
     
    180217                  zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
    181218                  !                                      ! uslp and vslp output in zwz and zww, resp. 
    182                   zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
    183                   zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
    184                   zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
    185                      &                   + zfi  * uslpml(ji,jj)                                                     & 
    186                      &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   & 
    187                      &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 
    188                   zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
    189                      &                   + zfj  * vslpml(ji,jj)                                                     & 
    190                      &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   & 
    191                      &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 
     219                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj  ,jk) ) 
     220                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) ) 
     221                  ! thickness of water column between surface and level k at u/v point 
     222                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              & 
     223                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) ) 
     224                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              & 
     225                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 
     226                  ! 
     227                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          & 
     228                     &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj) 
     229                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 
     230                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          & 
     231                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj)  
     232                  zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 
     233                   
     234                  
    192235!!gm  modif to suppress omlmask.... (as in Griffies case) 
    193236!                  !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0. 
     
    238281               DO ji = fs_2, fs_jpim1   ! vector opt. 
    239282                  uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
    240                      &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp 
     283                     &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   & 
     284                     &                            *   umask(ji,jj,jk-1) 
    241285                  vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
    242                      &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp 
     286                     &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   & 
     287                     &                            *   vmask(ji,jj,jk-1) 
    243288               END DO 
    244289            END DO 
     
    253298               DO ji = fs_2, fs_jpim1   ! vector opt. 
    254299                  !                                  !* Local vertical density gradient evaluated from N^2 
    255                   zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 
     300                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk) 
    256301                  !                                  !* Slopes at w point 
    257302                  !                                        ! i- & j-gradient of density at w-points 
     
    261306                     &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) *  e2t(ji,jj) 
    262307                  zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           & 
    263                      &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk) 
     308                     &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci 
    264309                  zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           & 
    265                      &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk) 
     310                     &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj 
    266311                  !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
    267312                  !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
     
    270315                  !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    271316                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
    272                   zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 
    273                   zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
    274                   zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
     317                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 
     318                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 
     319                     &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
     320                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 
     321                     &            + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    275322 
    276323!!gm  modif to suppress omlmask....  (as in Griffies operator) 
     
    325372                  zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
    326373                     &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
    327                   wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 
    328                   wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 
     374                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 
     375                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 
    329376               END DO 
    330377            END DO 
     
    377424         ! set the slope of diffusion to the slope of s-surfaces  
    378425         !      ( c a u t i o n : minus sign as fsdep has positive value )  
    379          DO jk = 1, jpk  
     426         DO jj = 2, jpjm1  
     427            DO ji = fs_2, fs_jpim1   ! vector opt.  
     428               uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) )  * umask(ji,jj,1)  
     429               vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) )  * vmask(ji,jj,1)  
     430               wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5  
     431               wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5  
     432            END DO  
     433         END DO  
     434 
     435         DO jk = 2, jpk  
    380436            DO jj = 2, jpjm1  
    381437               DO ji = fs_2, fs_jpim1   ! vector opt.  
    382438                  uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)  
    383439                  vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)  
    384                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5  
    385                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5  
     440                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 
     441                    &                              * wmask(ji,jj,jk) * 0.5  
     442                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 
     443                    &                              * wmask(ji,jj,jk) * 0.5  
    386444               END DO  
    387445            END DO  
     
    405463       
    406464      CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 
     465      CALL wrk_dealloc( jpi,jpj,     zhmlpu, zhmlpv) 
    407466      ! 
    408467      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp') 
     
    435494      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 
    436495      REAL(wp) ::   zdzrho_raw 
    437       REAL(wp) ::   zbeta0 
    438496      REAL(wp), POINTER, DIMENSION(:,:)     ::   z1_mlbw 
    439       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalbet 
    440497      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zdxrho , zdyrho, zdzrho     ! Horizontal and vertical density gradients 
    441498      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zti_mlb, ztj_mlb            ! for Griffies operator only 
     
    445502      ! 
    446503      CALL wrk_alloc( jpi,jpj, z1_mlbw ) 
    447       CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 
    448504      CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  ) 
    449505      CALL wrk_alloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  ) 
     
    452508      !  Some preliminary calculation  ! 
    453509      !--------------------------------! 
    454       ! 
    455       CALL eos_alpbet( tsb, zalbet, zbeta0 )  !==  before local thermal/haline expension ratio at T-points  ==! 
    456510      ! 
    457511      DO jl = 0, 1                            !==  unmasked before density i- j-, k-gradients  ==! 
     
    465519                  zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point 
    466520                  zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
    467                   zdxrho_raw = ( - zalbet(ji+ip,jj   ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 
    468                   zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 
    469                   zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
     521                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,jk,jp_tem) * zdit + rab_b(ji+ip,jj   ,jk,jp_sal) * zdis ) / e1u(ji,jj) 
     522                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji   ,jj+jp,jk,jp_sal) * zdjs ) / e2v(ji,jj) 
     523                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
    470524                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
    471525               END DO 
     
    473527         END DO 
    474528         ! 
    475          IF( ln_zps.and.l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom 
    476 # if defined key_vectopt_loop 
    477             DO jj = 1, 1 
    478                DO ji = 1, jpij-jpi            ! vector opt. (forced unrolling) 
    479 # else 
     529         IF( ln_zps .AND. l_grad_zps ) THEN     ! partial steps: correction of i- & j-grad on bottom 
    480530            DO jj = 1, jpjm1 
    481531               DO ji = 1, jpim1 
    482 # endif 
    483532                  iku  = mbku(ji,jj)          ;   ikv  = mbkv(ji,jj)             ! last ocean level (u- & v-points) 
    484533                  zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature 
    485534                  zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity 
    486                   zdxrho_raw = ( - zalbet(ji+ip,jj   ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 
    487                   zdyrho_raw = ( - zalbet(ji   ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 
     535                  zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) / e1u(ji,jj) 
     536                  zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) / e2v(ji,jj) 
    488537                  zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
    489538                  zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     
    505554                     zdks = 0._wp 
    506555                  ENDIF 
    507                   zdzrho_raw = ( - zalbet(ji   ,jj   ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 
    508                   zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )    ! force zdzrho >= repsln 
     556                  zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) * zdks ) / fse3w(ji,jj,jk+kp) 
     557                  zdzrho(ji,jj,jk,kp) = - MIN( - repsln, zdzrho_raw )    ! force zdzrho >= repsln 
    509558                 END DO 
    510559            END DO 
     
    650699      ! 
    651700      CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 
    652       CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 
    653701      CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho,              klstart = 0  ) 
    654702      CALL wrk_dealloc( jpi,jpj,  2,2, zti_mlb, ztj_mlb,        kkstart = 0, klstart = 0  ) 
     
    701749      !                                            !==   surface mixed layer mask   ! 
    702750      DO jk = 1, jpk                               ! =1 inside the mixed layer, =0 otherwise 
    703 # if defined key_vectopt_loop 
    704          DO jj = 1, 1 
    705             DO ji = 1, jpij                        ! vector opt. (forced unrolling) 
    706 # else 
    707751         DO jj = 1, jpj 
    708752            DO ji = 1, jpi 
    709 # endif 
    710753               ik = nmln(ji,jj) - 1 
    711                IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp 
    712                ELSE                  ;   omlmask(ji,jj,jk) = 0._wp 
     754               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN 
     755                  omlmask(ji,jj,jk) = 1._wp 
     756               ELSE 
     757                  omlmask(ji,jj,jk) = 0._wp 
    713758               ENDIF 
    714759            END DO 
     
    727772      !----------------------------------------------------------------------- 
    728773      ! 
    729 # if defined key_vectopt_loop 
    730       DO jj = 1, 1 
    731          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    732 # else 
    733774      DO jj = 2, jpjm1 
    734775         DO ji = 2, jpim1 
    735 # endif 
    736776            !                        !==   Slope at u- & v-points just below the Mixed Layer   ==! 
    737777            ! 
    738778            !                        !- vertical density gradient for u- and v-slopes (from dzr at T-point) 
    739             iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
    740             ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
     779            iku = MIN(  MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
     780            ikv = MIN(  MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
    741781            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) ) 
    742782            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) ) 
     
    772812            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  ) 
    773813            !                        !- i- & j-slope at w-points (wslpiml, wslpjml) 
    774             wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 
    775             wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 
     814            wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik) 
     815            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik) 
    776816         END DO 
    777817      END DO 
     
    825865         wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp 
    826866 
    827          IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN 
     867         IF(ln_sco .AND.  (ln_traldf_hor .OR. ln_dynldf_hor )) THEN 
    828868            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 
    829869 
Note: See TracChangeset for help on using the changeset viewer.