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 5189 for branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 – NEMO

Ignore:
Timestamp:
2015-03-31T19:58:23+02:00 (9 years ago)
Author:
mathiot
Message:

ISF cleaning branch: simplification and bug correction in hpg_isf, zps_hde_isf, mixed layer definition, slope, diffusion, vertical advection and top friction

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5120 r5189  
    196196   END SUBROUTINE zps_hde 
    197197   ! 
    198    SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    199       &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    200       &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
     198   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  & 
     199      &                          prd, pgru, pgrv, pgrui, pgrvi ) 
    201200      !!---------------------------------------------------------------------- 
    202201      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    252251      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    253252      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    254       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom) 
    255       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    256       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    257253      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top) 
    258       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top) 
    259       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top) 
    260       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
    261254      ! 
    262255      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    … …  
    269262      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    270263      ! 
    271       pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
     264      pgtu (:,:,:)=0.0_wp ; pgtv(:,:,:) =0.0_wp ; 
    272265      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    273       zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    274       zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     266      zti  (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
     267      zhi  (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
    275268      ! 
    276269      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    … …  
    280273               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    281274               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     275               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
     276               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     277               ! 
     278               ! i- direction 
     279               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     280                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     281                  ! interpolated values of tracers 
     282                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     283                  ! gradient of  tracers 
     284                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     285               ELSE                           ! case 2 
     286                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     287                  ! interpolated values of tracers 
     288                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     289                  ! gradient of tracers 
     290                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     291               ENDIF 
     292               ! 
     293               ! j- direction 
     294               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     295                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     296                  ! interpolated values of tracers 
     297                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     298                  ! gradient of tracers 
     299                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     300               ELSE                           ! case 2 
     301                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     302                  ! interpolated values of tracers 
     303                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     304                  ! gradient of tracers 
     305                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     306               ENDIF 
     307            END DO 
     308         END DO 
     309         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     310         ! 
     311      END DO 
     312 
     313      ! horizontal derivative of density anomalies (rd) 
     314      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     315         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     316         ! 
     317         DO jj = 1, jpjm1 
     318            DO ji = 1, jpim1 
     319               iku = mbku(ji,jj) 
     320               ikv = mbkv(ji,jj) 
     321               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
     322               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     323 
     324               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
     325               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
     326               ENDIF 
     327               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
     328               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
     329               ENDIF 
     330 
     331            END DO 
     332         END DO 
     333 
     334         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     335         ! step and store it in  zri, zrj for each  case 
     336         CALL eos( zti, zhi, zri ) 
     337         CALL eos( ztj, zhj, zrj ) 
     338 
     339         ! Gradient of density at the last level  
     340         DO jj = 1, jpjm1 
     341            DO ji = 1, jpim1 
     342               iku = mbku(ji,jj) 
     343               ikv = mbkv(ji,jj) 
     344               ze3wu = gdept_0(ji+1,jj,iku) - gdept_0(ji,jj,iku) 
     345               ze3wv = gdept_0(ji,jj+1,ikv) - gdept_0(ji,jj,ikv) 
     346 
     347               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     348               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     349               ENDIF 
     350               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     351               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     352               ENDIF 
     353            END DO 
     354         END DO 
     355         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     356         ! 
     357      END IF 
     358         ! (ISH)  compute grui and gruvi 
     359      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
     360         DO jj = 1, jpjm1 
     361            DO ji = 1, jpim1 
     362               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
     363               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
     364               ! 
    282365               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    283366               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    284367               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    285368               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    286                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    287                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    288                ! 
     369               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
     370               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
    289371               ! i- direction 
    290372               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    291                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    292                   ! interpolated values of tracers 
    293                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     373                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
     374                  ! interpolated values of tracers 
     375                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
     376                  ! gradient of tracers 
     377                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     378               ELSE                           ! case 2 
     379                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
     380                  ! interpolated values of tracers 
     381                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    294382                  ! gradient of  tracers 
    295                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    296                ELSE                           ! case 2 
    297                   zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    298                   ! interpolated values of tracers 
    299                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    300                   ! gradient of tracers 
    301                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     383                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    302384               ENDIF 
    303385               ! 
    304386               ! j- direction 
    305387               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    306                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    307                   ! interpolated values of tracers 
    308                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    309                   ! gradient of tracers 
    310                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    311                ELSE                           ! case 2 
    312                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    313                   ! interpolated values of tracers 
    314                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    315                   ! gradient of tracers 
    316                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    317                ENDIF 
    318             END DO 
    319          END DO 
    320          CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     388                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
     389                  ! interpolated values of tracers 
     390                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
     391                  ! gradient of tracers 
     392                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     393               ELSE                           ! case 2 
     394                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
     395                  ! interpolated values of tracers 
     396                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
     397                  ! gradient of tracers 
     398                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     399               ENDIF 
     400            END DO!! 
     401         END DO!! 
     402         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    321403         ! 
    322404      END DO 
    … …  
    324406      ! horizontal derivative of density anomalies (rd) 
    325407      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    326          pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
    327          pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ; 
    328          pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ; 
    329          pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ; 
    330          DO jj = 1, jpjm1 
    331             DO ji = 1, jpim1 
    332                iku = mbku(ji,jj) 
    333                ikv = mbkv(ji,jj) 
    334                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    335                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    336  
    337                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1 
    338                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2 
    339                ENDIF 
    340                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1 
    341                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2 
    342                ENDIF 
    343             END DO 
    344          END DO 
    345           
     408         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
     409         DO jj = 1, jpjm1 
     410            DO ji = 1, jpim1 
     411               iku = miku(ji,jj) 
     412               ikv = mikv(ji,jj) 
     413               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
     414               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
     415 
     416               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
     417               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
     418               ENDIF 
     419 
     420               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
     421               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
     422               ENDIF 
     423 
     424            END DO 
     425         END DO 
     426 
    346427         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    347428         ! step and store it in  zri, zrj for each  case 
    … …  
    352433         DO jj = 1, jpjm1 
    353434            DO ji = 1, jpim1 
    354                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    355                ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    356                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    357                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    358                IF( ze3wu >= 0._wp ) THEN  
    359                   pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 
    360                   pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    361                   pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1  
    362                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    363                                 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) & 
    364                                    - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    365                ELSE   
    366                   pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 
    367                   pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    368                   pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2 
    369                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    370                                 * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 
    371                                    -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    372                ENDIF 
    373                IF( ze3wv >= 0._wp ) THEN 
    374                   pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv)  
    375                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    376                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1 
    377                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    378                                 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    379                                    - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    380                ELSE  
    381                   pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 
    382                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    383                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2 
    384                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    385                                 * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    386                                    -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    387                ENDIF 
    388             END DO 
    389          END DO 
    390          CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions 
    391          CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions 
    392          CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions 
    393          CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions 
    394          ! 
    395       END IF 
    396          ! (ISH)  compute grui and gruvi 
    397       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
    398          DO jj = 1, jpjm1 
    399             DO ji = 1, jpim1 
    400                iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
    401                ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
    402                ! 
    403                ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    404                ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    405                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    406                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    407                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))  
    408                ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    409                ! i- direction 
    410                IF( ze3wu >= 0._wp ) THEN      ! case 1 
    411                   zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
    412                   ! interpolated values of tracers 
    413                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    414                   ! gradient of tracers 
    415                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    416                ELSE                           ! case 2 
    417                   zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    418                   ! interpolated values of tracers 
    419                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    420                   ! gradient of  tracers 
    421                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    422                ENDIF 
    423                ! 
    424                ! j- direction 
    425                IF( ze3wv >= 0._wp ) THEN      ! case 1 
    426                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
    427                   ! interpolated values of tracers 
    428                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    429                   ! gradient of tracers 
    430                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    431                ELSE                           ! case 2 
    432                   zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    433                   ! interpolated values of tracers 
    434                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    435                   ! gradient of tracers 
    436                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    437                ENDIF 
    438             END DO!! 
    439          END DO!! 
    440          CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    441          ! 
    442       END DO 
    443  
    444       ! horizontal derivative of density anomalies (rd) 
    445       IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    446          pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
    447          pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
    448          pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
    449          pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    450  
    451          DO jj = 1, jpjm1 
    452             DO ji = 1, jpim1 
    453                iku = miku(ji,jj) 
    454                ikv = mikv(ji,jj) 
    455                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    456                ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    457  
    458                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    459                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    460                ENDIF 
    461                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
    462                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
    463                ENDIF 
    464             END DO 
    465          END DO 
    466  
    467          ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    468          ! step and store it in  zri, zrj for each  case 
    469          CALL eos( zti, zhi, zri )   
    470          CALL eos( ztj, zhj, zrj ) 
    471  
    472          ! Gradient of density at the last level  
    473          DO jj = 1, jpjm1 
    474             DO ji = 1, jpim1 
    475435               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
    476436               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 
    477                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    478                ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    479                IF( ze3wu >= 0._wp ) THEN 
    480                  pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
    481                  pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
    482                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    483                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    484                                 * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    485                                    - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    486                ELSE 
    487                  pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    488                  pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
    489                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    490                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    491                                 * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    492                                    -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    493                ENDIF 
    494                IF( ze3wv >= 0._wp ) THEN 
    495                  pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
    496                  pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
    497                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    498                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    499                                 * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    500                                    - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    501                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    502                ELSE 
    503                  pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
    504                  pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
    505                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    506                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    507                                 * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    508                                    -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    509                ENDIF 
     437               ze3wu  =  gdept_0(ji,jj,iku) - gdept_0(ji+1,jj,iku) 
     438               ze3wv  =  gdept_0(ji,jj,ikv) - gdept_0(ji,jj+1,ikv)  
     439 
     440               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) )          ! i: 1 
     441               ELSE                      ; pgrui(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) )      ! i: 2 
     442               ENDIF 
     443 
     444               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1 
     445               ELSE                      ; pgrvi(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2 
     446               ENDIF 
     447 
    510448            END DO 
    511449         END DO 
    512450         CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
    513          CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions 
    514          CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions 
    515          CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions 
    516451         ! 
    517452      END IF   
Note: See TracChangeset for help on using the changeset viewer.