Changeset 2840


Ignore:
Timestamp:
2011-09-15T14:38:15+02:00 (9 years ago)
Author:
agn
Message:

branches/2011/dev_r2782_NOCS_Griffies ticket #838 bugfixes + preliminary support for s-coords

Location:
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO
Files:
6 edited

Legend:

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

    r2772 r2840  
    4646   !                                                                !! Griffies operator 
    4747   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)     ::   wslp2                !: wslp**2 from Griffies quarter cells 
    48    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials  
     48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials 
    4949   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate 
    5050 
     
    5858 
    5959   ! Workspace arrays for ldf_slp_grif. These could be replaced by several 3D and 2D workspace 
    60    ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace  
     60   ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace 
    6161   ! arrays can't be used here because of the zero-indexing of some of the ranks. ARPDBG. 
    6262   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   zdzrho , zdyrho, zdxrho     ! Horizontal and vertical density gradients 
     
    9393      !!---------------------------------------------------------------------- 
    9494      !!                 ***  ROUTINE ldf_slp  *** 
    95       !!  
     95      !! 
    9696      !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal 
    9797      !!              surfaces referenced locally) (ln_traldf_iso=T). 
    98       !!  
    99       !! ** Method  :   The slope in the i-direction is computed at U- and  
    100       !!      W-points (uslp, wslpi) and the slope in the j-direction is  
     98      !! 
     99      !! ** Method  :   The slope in the i-direction is computed at U- and 
     100      !!      W-points (uslp, wslpi) and the slope in the j-direction is 
    101101      !!      computed at V- and W-points (vslp, wslpj). 
    102102      !!      They are bounded by 1/100 over the whole ocean, and within the 
     
    112112      !!      bottom slope (ln_sco=T) at level jpk in inildf] 
    113113      !! 
    114       !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes  
     114      !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and  j-slopes 
    115115      !!               of now neutral surfaces at u-, w- and v- w-points, resp. 
    116116      !!---------------------------------------------------------------------- 
     
    127127      INTEGER  ::   ii0, ii1, iku   ! temporary integer 
    128128      INTEGER  ::   ij0, ij1, ikv   ! temporary integer 
    129       REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16    ! local scalars 
     129      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars 
    130130      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      - 
    131131      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
     
    148148         DO jj = 1, jpjm1 
    149149            DO ji = 1, fs_jpim1   ! vector opt. 
    150                zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) )  
    151                zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) )  
     150               zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj  ,jk) - prd(ji,jj,jk) ) 
     151               zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji  ,jj+1,jk) - prd(ji,jj,jk) ) 
    152152            END DO 
    153153         END DO 
    154154      END DO 
    155155      IF( ln_zps ) THEN                           ! partial steps correction at the bottom ocean level 
    156 # if defined key_vectopt_loop   
     156# if defined key_vectopt_loop 
    157157         DO jj = 1, 1 
    158158            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     
    161161            DO ji = 1, jpim1 
    162162# endif 
    163                zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)  
    164                zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)                
     163               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
     164               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
    165165            END DO 
    166166         END DO 
     
    181181      CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr )        ! output: uslpml, vslpml, wslpiml, wslpjml 
    182182 
    183        
     183 
    184184      ! I.  slopes at u and v point      | uslp = d/di( prd ) / d/dz( prd ) 
    185185      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    186       !                
     186      ! 
    187187      DO jk = 2, jpkm1                            !* Slopes at u and v points 
    188188         DO jj = 2, jpjm1 
     
    225225      DO jk = 2, jpkm1 
    226226         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only 
    227             DO ji = 2, jpim1   
     227            DO ji = 2, jpim1 
    228228               uslp(ji,jj,jk) = z1_16 * (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)      & 
    229229                  &                       +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)      & 
     
    266266      ! II.  slopes at w point           | wslpi = mij( d/di( prd ) / d/dz( prd ) 
    267267      ! ===========================      | wslpj = mij( d/dj( prd ) / d/dz( prd ) 
    268       !                
     268      ! 
    269269      DO jk = 2, jpkm1 
    270270         DO jj = 2, jpjm1 
     
    308308         DO jj = 2, jpjm1, MAX(1, jpj-3)                        ! rows jj=2 and =jpjm1 only 
    309309            DO ji = 2, jpim1 
     310               zcofw = tmask(ji,jj,jk) * z1_16 
    310311               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
    311                   &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
    312                   &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
    313                   &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
    314                   &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
     312                    &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
     313                    &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
     314                    &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
     315                    &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw 
    315316 
    316317               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
    317                   &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
    318                   &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
    319                   &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
    320                   &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
    321             END DO 
    322          END DO   
     318                    &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
     319                    &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
     320                    &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
     321                    &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw 
     322            END DO 
     323         END DO 
    323324         DO jj = 3, jpj-2                               ! other rows 
    324325            DO ji = fs_2, fs_jpim1   ! vector opt. 
     326               zcofw = tmask(ji,jj,jk) * z1_16 
    325327               wslpi(ji,jj,jk) = (        zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk)     & 
    326                   &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
    327                   &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
    328                   &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
    329                   &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
     328                    &                +      zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk)     & 
     329                    &                + 2.*( zwz(ji  ,jj-1,jk) + zwz(ji-1,jj  ,jk)     & 
     330                    &                +      zwz(ji+1,jj  ,jk) + zwz(ji  ,jj+1,jk) )   & 
     331                    &                + 4.*  zwz(ji  ,jj  ,jk)                         ) * zcofw 
    330332 
    331333               wslpj(ji,jj,jk) = (        zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk)     & 
    332                   &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
    333                   &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
    334                   &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
    335                   &                + 4.*  zww(ji  ,jj  ,jk)                         ) * z1_16 * tmask(ji,jj,jk) 
     334                    &                +      zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk)     & 
     335                    &                + 2.*( zww(ji  ,jj-1,jk) + zww(ji-1,jj  ,jk)     & 
     336                    &                +      zww(ji+1,jj  ,jk) + zww(ji  ,jj+1,jk) )   & 
     337                    &                + 4.*  zww(ji  ,jj  ,jk)                         ) * zcofw 
    336338            END DO 
    337339         END DO 
     
    346348         END DO 
    347349      END DO 
    348           
    349       ! III.  Specific grid points      
    350       ! ===========================  
    351       !                
     350 
     351      ! III.  Specific grid points 
     352      ! =========================== 
     353      ! 
    352354      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN     !  ORCA_R4 configuration: horizontal diffusion in specific area 
    353355         !                                                    ! Gibraltar Strait 
     
    368370      ENDIF 
    369371 
    370       ! IV. Lateral boundary conditions  
     372 
     373      ! IV. Lateral boundary conditions 
    371374      ! =============================== 
    372375      CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
     
    382385      ! 
    383386   END SUBROUTINE ldf_slp 
    384     
     387 
    385388 
    386389   SUBROUTINE ldf_slp_grif ( kt ) 
     
    390393      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope 
    391394      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 
    392       !!      at W-points using the Griffies quarter-cells.   
    393       !! 
    394       !! ** Method  :   calculates alpha and beta at T-points  
     395      !!      at W-points using the Griffies quarter-cells. 
     396      !! 
     397      !! ** Method  :   calculates alpha and beta at T-points 
    395398      !! 
    396399      !! ** Action : - triadi_g, triadj_g   T-pts i- and j-slope triads relative to geopot. (used for eiv) 
     
    399402      !!---------------------------------------------------------------------- 
    400403      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    401       USE oce     , ONLY:   zdit    => ua       , zdis   => va         ! (ua,va) used as workspace 
    402       USE oce     , ONLY:   zdjt    => ta       , zdjs   => sa         ! (ta,sa) used as workspace 
    403       USE wrk_nemo, ONLY:   zdkt    => wrk_3d_2 , zdks   => wrk_3d_3   ! 3D workspace 
    404       USE wrk_nemo, ONLY:   zalpha  => wrk_3d_4 , zbeta => wrk_3d_5    ! alpha, beta at T points, at depth fsgdept 
     404      USE oce     , ONLY:   zalbet  => ua   ! use ua as workspace 
    405405      USE wrk_nemo, ONLY:   z1_mlbw => wrk_2d_1 
    406       ! 
    407       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    408       ! 
     406      !! 
     407      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
     408      !! 
    409409      INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices 
    410       INTEGER  ::   iku, ikv                                  ! local integer 
    411       REAL(wp) ::   zfacti, zfactj, zatempw,zatempu,zatempv   ! local scalars 
    412       REAL(wp) ::   zbu, zbv, zbti, zbtj                      !   -      - 
     410      INTEGER  ::   iku, ikv                    ! local integer 
     411      REAL(wp) ::   zfacti, zfactj              ! local scalars 
     412      REAL(wp) ::   zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 
    413413      REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 
    414414      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 
    415415      REAL(wp) ::   zdzrho_raw 
     416      REAL(wp) ::   zbeta0 
    416417      !!---------------------------------------------------------------------- 
    417418 
     
    424425      !--------------------------------! 
    425426      ! 
    426       CALL eos_alpbet( tsb, zalpha, zbeta )     !==  before thermal and haline expension coeff. at T-points  ==! 
    427       ! 
    428       DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==! 
    429          DO jj = 1, jpjm1 
    430             DO ji = 1, fs_jpim1   ! vector opt. 
    431                zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk)   ! i-gradient of T and S at jj 
    432                zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 
    433                zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk)   ! j-gradient of T and S at jj 
    434                zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 
    435             END DO 
    436          END DO 
    437       END DO 
    438       IF( ln_zps ) THEN                               ! partial steps: correction at the last level 
     427      CALL eos_alpbet( tsb, zalbet, zbeta0 )         !==  before local thermal/haline expension ratio at T-points  ==! 
     428      ! 
     429      DO jl = 0, 1                           !==  unmasked before density i- j-, k-gradients  ==! 
     430         ! 
     431         ip = jl   ;   jp = jl        ! guaranteed nonzero gradients ( absolute value larger than repsln) 
     432         DO jk = 1, jpkm1                          ! done each pair of triad 
     433            DO jj = 1, jpjm1                       ! NB: not masked ==>  a minimum value is set 
     434               DO ji = 1, fs_jpim1   ! vector opt. 
     435                  zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! i-gradient of T & S at u-point 
     436                  zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     437                  zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )    ! j-gradient of T & S at v-point 
     438                  zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 
     439                  zdxrho_raw = ( - zalbet(ji+ip,jj   ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 
     440                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 
     441                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
     442                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     443               END DO 
     444            END DO 
     445         END DO 
     446         ! 
     447         IF( ln_zps.and.ln_grad_zps ) THEN                       ! partial steps: correction of i- & j-grad on bottom 
    439448# if defined key_vectopt_loop 
    440          DO jj = 1, 1 
    441             DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     449            DO jj = 1, 1 
     450               DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    442451# else 
    443          DO jj = 1, jpjm1 
    444             DO ji = 1, jpim1 
     452            DO jj = 1, jpjm1 
     453               DO ji = 1, jpim1 
    445454# endif 
    446                zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem)                           ! i-gradient of T and S 
    447                zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal) 
    448                zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem)                           ! j-gradient of T and S 
    449                zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal) 
    450             END DO 
    451          END DO 
    452       ENDIF 
    453       ! 
    454       zdkt(:,:,1) = 0._wp                    !==  before vertical T & S gradient at w-level  ==! 
    455       zdks(:,:,1) = 0._wp 
    456       DO jk = 2, jpk 
    457          zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 
    458          zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 
    459       END DO 
    460       ! 
    461       ! 
    462       DO jl = 0, 1                           !==  density i-, j-, and k-gradients  ==! 
    463          ip = jl   ;   jp = jl         ! guaranteed nonzero gradients ( absolute value larger than repsln) 
     455                  iku  = mbku(ji,jj)          ;   ikv  = mbkv(ji,jj)             ! last ocean level (u- & v-points) 
     456                  zdit = gtsu(ji,jj,jp_tem)   ;   zdjt = gtsv(ji,jj,jp_tem)      ! i- & j-gradient of Temperature 
     457                  zdis = gtsu(ji,jj,jp_sal)   ;   zdjs = gtsv(ji,jj,jp_sal)      ! i- & j-gradient of Salinity 
     458                  zdxrho_raw = ( - zalbet(ji+ip,jj   ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 
     459                  zdyrho_raw = ( - zalbet(ji   ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 
     460                  zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
     461                  zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     462               END DO 
     463            END DO 
     464         ENDIF 
     465         ! 
     466      END DO 
     467 
     468      DO kp = 0, 1                           !==  unmasked before density i- j-, k-gradients  ==! 
    464469         DO jk = 1, jpkm1                          ! done each pair of triad 
    465             DO jj = 1, jpjm1                       ! NB: not masked due to the minimum value set 
    466                DO ji = 1, fs_jpim1   ! vector opt.  
    467                   zdxrho_raw = ( zalpha(ji+ip,jj   ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj   ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj) 
    468                   zdyrho_raw = ( zalpha(ji   ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji   ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj) 
    469                   zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )    ! keep the sign 
    470                   zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX(   repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
    471                END DO 
    472             END DO 
    473          END DO 
    474       END DO 
    475      DO kp = 0, 1                           !==  density i-, j-, and k-gradients  ==! 
    476          DO jk = 1, jpkm1                          ! done each pair of triad 
    477             DO jj = 1, jpj                       ! NB: not masked due to the minimum value set 
    478                DO ji = 1, jpi   ! vector opt.  
    479                   zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) )   & 
    480                      &       / fse3w(ji,jj,jk+kp) 
    481                   zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )                    ! force zdzrho >= repsln 
    482                END DO 
     470            DO jj = 1, jpj                       ! NB: not masked ==>  a minimum value is set 
     471               DO ji = 1, jpi  ! vector opt. 
     472                  IF( jk+kp > 1 ) THEN                                       ! k-gradient of T & S a jk+kp 
     473                     zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 
     474                     zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 
     475                  ELSE 
     476                     zdkt = 0._wp                                             ! 1st level gradient set to zero 
     477                     zdks = 0._wp 
     478                  ENDIF 
     479                  zdzrho_raw = ( - zalbet(ji   ,jj   ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 
     480                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )    ! force zdzrho >= repsln 
     481                 END DO 
    483482            END DO 
    484483         END DO 
     
    494493      !                                      !==  intialisations to zero  ==! 
    495494      ! 
    496       wslp2  (:,:,:)     = 0._wp                                           ! wslp2 will be cumulated 3D field set to zero 
    497       triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp      ! set surface and bottom slope to zero 
     495      wslp2  (:,:,:)     = 0._wp                                        ! wslp2 will be cumulated 3D field set to zero 
     496      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero 
    498497      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp 
    499498!!gm _iso set to zero missing 
    500       triadi (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero 
    501       triadj (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp 
    502        
     499      triadi  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero 
     500      triadj  (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp 
     501 
    503502      !-------------------------------------! 
    504503      !  Triads just below the Mixed Layer  ! 
     
    506505      ! 
    507506      DO jl = 0, 1               ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 
    508          DO kp = 0, 1            ! with only the slope-max limit   and   MASKED  
     507         DO kp = 0, 1            ! with only the slope-max limit   and   MASKED 
    509508            DO jj = 1, jpjm1 
    510509               DO ji = 1, fs_jpim1 
    511510                  ip = jl   ;   jp = jl 
    512511                  jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1         ! ML level+1 (MIN in case ML depth is the ocean depth) 
     512                  ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 
    513513                  zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      & 
    514                      &      + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk) 
     514                     &      - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk) 
    515515                  jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 
    516516                  ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      & 
    517                      &      + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk) 
     517                     &      - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk) 
    518518                  zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 
    519519                  ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 
     
    535535                     ! 
    536536                     ! Calculate slope relative to geopotentials used for GM skew fluxes 
    537                      ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth) 
     537                     ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 
    538538                     ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 
    539539                     ! masked by umask taken at the level of dz(rho) 
     
    544544                     ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp) 
    545545                     zti_coord = ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
    546                      ztj_coord = ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
    547                   ! unmasked 
    548                      zti_g_raw = zti_raw + zti_coord      ! ref to geopot surfaces 
    549                      ztj_g_raw = ztj_raw + ztj_coord 
     546                     ztj_coord = ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)                  ! unmasked 
     547                     zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces 
     548                     ztj_g_raw = ztj_raw - ztj_coord 
    550549                     zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 
    551550                     ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 
     
    562561                     ztj_g_lim =           zfactj   * ztj_g_lim                       & 
    563562                        &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   & 
    564                         &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp)                   ! masked 
    565                      ! 
    566                      triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) 
     563                        &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) 
     564                     ! 
     565                     triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp)                     ! masked 
    567566                     triadj_g(ji   ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp) 
    568567                     ! 
     
    574573                     ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 
    575574                     ! 
    576                      zti_lim  = zti_g_lim - zti_coord    ! remove the coordinate slope  ==> relative to coordinate surfaces 
    577                      ztj_lim  = ztj_g_lim - ztj_coord                                  
    578                      zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp)      ! square of limited slopes            ! masked <<== 
    579                      ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp) 
    580                      ! 
     575                     zti_lim  = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp)    ! remove coordinate slope => relative to coordinate surfaces 
     576                     ztj_lim  = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 
     577                     zti_lim2 = zti_lim * zti_lim      ! square of limited slopes            ! masked <<== 
     578                     ztj_lim2 = ztj_lim * ztj_lim 
     579                     ! 
     580                     IF( ln_triad_iso ) THEN 
     581                        zti_raw = zti_lim2 / zti_raw 
     582                        ztj_raw = ztj_lim2 / ztj_raw 
     583                        zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 
     584                        ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 
     585                        zti_lim =           zfacti   * zti_lim                       & 
     586                        &      + ( 1._wp - zfacti ) * zti_raw 
     587                        ztj_lim =           zfactj   * ztj_lim                       & 
     588                        &      + ( 1._wp - zfactj ) * ztj_raw 
     589                     ENDIF 
     590                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim 
     591                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim 
     592                    ! 
    581593                     zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   ) 
    582594                     zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   ) 
     
    584596                     zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 
    585597                     ! 
    586                      triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim2 / zti_raw                                          ! masked 
    587                      triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw 
    588                      ! 
    589598!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked 
    590                      wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2             ! masked 
     599! agn may need to change to using ztj_g_lim**2, as wslp2 just used for eddy growth rate, needs *real* slope 
     600                     wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2      ! masked 
    591601                     wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 
    592602                  END DO 
     
    597607      ! 
    598608      wslp2(:,:,1) = 0._wp                ! force the surface wslp to zero 
    599        
     609 
    600610      CALL lbc_lnk( wslp2, 'W', 1. )      ! lateral boundary confition on wslp2 only   ==>>> gm : necessary ? to be checked 
    601611      ! 
     
    610620      !!                  ***  ROUTINE ldf_slp_mxl  *** 
    611621      !! 
    612       !! ** Purpose :   Compute the slopes of iso-neutral surface just below  
     622      !! ** Purpose :   Compute the slopes of iso-neutral surface just below 
    613623      !!              the mixed layer. 
    614624      !! 
     
    619629      !! 
    620630      !! ** Action  :   uslpml, wslpiml :  i- &  j-slopes of neutral surfaces 
    621       !!                vslpml, wslpjml    just below the mixed layer  
     631      !!                vslpml, wslpjml    just below the mixed layer 
    622632      !!                omlmask         :  mixed layer mask 
    623633      !!---------------------------------------------------------------------- 
     
    683693            !                          !- vertical density gradient for u- and v-slopes (from dzr at T-point) 
    684694            iku = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1  )   ! ML (MAX of T-pts, bound by jpkm1) 
    685             ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   !  
     695            ikv = MIN(  MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1  )   ! 
    686696            zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj  ,iku) ) 
    687697            zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji  ,jj+1,ikv) ) 
     
    705715            !                          !- horizontal density i- & j-gradient at w-points 
    706716            zci = MAX(   umask(ji-1,jj,ik  ) + umask(ji,jj,ik  )           & 
    707                &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj)  
     717               &       + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps  ) * e1t(ji,jj) 
    708718            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           & 
    709719               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj) 
     
    734744      !! ** Purpose :   Initialization for the isopycnal slopes computation 
    735745      !! 
    736       !! ** Method  :   read the nammbf namelist and check the parameter  
    737       !!              values called by tra_dmp at the first timestep (nit000) 
     746      !! ** Method  :   read the nammbf namelist and check the parameter 
     747      !!      values called by tra_dmp at the first timestep (nit000) 
    738748      !!---------------------------------------------------------------------- 
    739749      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    740750      INTEGER ::   ierr         ! local integer 
    741751      !!---------------------------------------------------------------------- 
    742        
    743       IF(lwp) THEN     
     752 
     753      IF(lwp) THEN 
    744754         WRITE(numout,*) 
    745755         WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 
    746756         WRITE(numout,*) '~~~~~~~~~~~~' 
    747757      ENDIF 
    748        
     758 
    749759      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes 
    750760         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) 
     
    755765         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 
    756766         ! 
    757          IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco )   & 
    758             CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' ) 
     767         ! IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco )   & 
     768         !    CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' ) 
    759769         ! 
    760770      ELSE                             ! Madec operator : slopes at u-, v-, and w-points 
     
    804814CONTAINS 
    805815   SUBROUTINE ldf_slp( kt, prd, pn2 )        ! Dummy routine 
    806       INTEGER, INTENT(in) :: kt  
     816      INTEGER, INTENT(in) :: kt 
    807817      REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 
    808818      WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) 
  • branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90

    r2715 r2840  
    6767         &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   & 
    6868         &                 ln_traldf_grif , ln_traldf_gdia,                   & 
     69         &                 ln_triad_iso , ln_botmix, ln_grad_zps,    & 
    6970         &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       & 
    7071         &                 rn_slpmax 
     
    9495         WRITE(numout,*) '      maximum isoppycnal slope      rn_slpmax       = ', rn_slpmax 
    9596         WRITE(numout,*) '   + griffies operator internal controls not set via the namelist (experimental): ' 
    96          WRITE(numout,*) '      calculate triads twice        l_triad_iso     = ', l_triad_iso 
    97          WRITE(numout,*) '      no Shapiro filter             l_no_smooth     = ', l_no_smooth 
     97         WRITE(numout,*) '      calculate triads twice        ln_triad_iso     = ', ln_triad_iso 
     98         WRITE(numout,*) '      special Tgradts w ptl steps   ln_grad_zps     = ', ln_grad_zps 
     99         WRITE(numout,*) '      GM -->lat mixing on bottom    ln_botmix       = ', ln_botmix 
    98100         WRITE(numout,*) 
    99101      ENDIF 
  • branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90

    r2715 r2840  
    3232 
    3333   REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names 
    34    LOGICAL , PUBLIC ::   l_triad_iso     = .FALSE.   !: calculate triads twice 
    35    LOGICAL , PUBLIC ::   l_no_smooth     = .FALSE.   !: no Shapiro smoothing 
     34   LOGICAL , PUBLIC ::   ln_triad_iso     = .FALSE.   !: calculate triads twice 
     35   LOGICAL , PUBLIC ::   ln_grad_zps     = .FALSE.   !: special treatment for Horz Tgradients w partial steps 
     36   LOGICAL , PUBLIC ::   ln_botmix     = .FALSE.   !:  mixing on bottom 
    3637 
    3738#if defined key_traldf_c3d 
  • branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r2715 r2840  
    33   !!                       ***  MODULE  eosbn2  *** 
    44   !! Ocean diagnostic variable : equation of state - in situ and potential density 
    5    !!                                               - Brunt-Vaisala frequency  
     5   !!                                               - Brunt-Vaisala frequency 
    66   !!============================================================================== 
    77   !! History :  OPA  ! 1989-03  (O. Marti)  Original code 
     
    2727   !!   eos_insitu_2d  : Compute the in situ density for 2d fields 
    2828   !!   eos_bn2        : Compute the Brunt-Vaisala frequency 
    29    !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff. 
     29   !!   eos_alpbet     : calculates the in situ thermal/haline expansion ratio 
    3030   !!   tfreez         : Compute the surface freezing temperature 
    3131   !!   eos_init       : set eos parameters (namelist) 
     
    4141   PRIVATE 
    4242 
    43    !                   !! * Interface  
     43   !                   !! * Interface 
    4444   INTERFACE eos 
    4545      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 
    46    END INTERFACE  
     46   END INTERFACE 
    4747   INTERFACE bn2 
    4848      MODULE PROCEDURE eos_bn2 
    49    END INTERFACE  
     49   END INTERFACE 
    5050 
    5151   PUBLIC   eos        ! called by step, istate, tranpc and zpsgrd modules 
     
    6161 
    6262   REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio 
    63     
     63 
    6464   !! * Substitutions 
    6565#  include "domzgr_substitute.h90" 
     
    7575      !!---------------------------------------------------------------------- 
    7676      !!                   ***  ROUTINE eos_insitu  *** 
    77       !!  
    78       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from  
     77      !! 
     78      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    7979      !!       potential temperature and salinity using an equation of state 
    8080      !!       defined through the namelist parameter nn_eos. 
     
    134134!CDIR NOVERRCHK 
    135135         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    136          !   
     136         ! 
    137137         DO jk = 1, jpkm1 
    138138            DO jj = 1, jpj 
     
    199199      !!---------------------------------------------------------------------- 
    200200      !!                  ***  ROUTINE eos_insitu_pot  *** 
    201       !!            
     201      !! 
    202202      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) and the 
    203203      !!      potential volumic mass (Kg/m3) from potential temperature and 
    204       !!      salinity fields using an equation of state defined through the  
     204      !!      salinity fields using an equation of state defined through the 
    205205      !!     namelist parameter nn_eos. 
    206206      !! 
     
    230230      !!      nn_eos = 2 : linear equation of state function of temperature and 
    231231      !!               salinity 
    232       !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0  
     232      !!              prd(t,s) = ( rho(t,s) - rau0 ) / rau0 
    233233      !!                       = rn_beta * s - rn_alpha * tn - 1. 
    234234      !!              rhop(t,s)  = rho(t,s) 
     
    265265!CDIR NOVERRCHK 
    266266         zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 
    267          !   
     267         ! 
    268268         DO jk = 1, jpkm1 
    269269            DO jj = 1, jpj 
     
    336336      !!                  ***  ROUTINE eos_insitu_2d  *** 
    337337      !! 
    338       !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from  
     338      !! ** Purpose :   Compute the in situ density (ratio rho/rau0) from 
    339339      !!      potential temperature and salinity using an equation of state 
    340340      !!      defined through the namelist parameter nn_eos. * 2D field case 
     
    374374      !                                                           ! 2 : salinity               [psu] 
    375375      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                  [m] 
    376       REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density  
     376      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density 
    377377      !! 
    378378      INTEGER  ::   ji, jj                    ! dummy loop indices 
     
    449449         DO jj = 1, jpjm1 
    450450            DO ji = 1, fs_jpim1   ! vector opt. 
    451                prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)  
     451               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
    452452            END DO 
    453453         END DO 
     
    468468      !! ** Purpose :   Compute the local Brunt-Vaisala frequency at the time- 
    469469      !!      step of the input arguments 
    470       !!       
     470      !! 
    471471      !! ** Method : 
    472472      !!       * nn_eos = 0  : UNESCO sea water properties 
     
    482482      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
    483483      !!      The use of potential density to compute N^2 introduces e r r o r 
    484       !!      in the sign of N^2 at great depths. We recommand the use of  
     484      !!      in the sign of N^2 at great depths. We recommand the use of 
    485485      !!      nn_eos = 0, except for academical studies. 
    486486      !!        Macro-tasked on horizontal slab (jk-loop) 
     
    497497      !! 
    498498      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    499       REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars  
     499      REAL(wp) ::   zgde3w, zt, zs, zh, zalbet, zbeta   ! local scalars 
    500500#if defined key_zdfddm 
    501501      REAL(wp) ::   zds   ! local scalars 
     
    504504 
    505505      ! pn2 : interior points only (2=< jk =< jpkm1 ) 
    506       ! --------------------------  
     506      ! -------------------------- 
    507507      ! 
    508508      SELECT CASE( nn_eos ) 
     
    542542                     &                                - 0.121555e-07_wp ) * zh 
    543543                     ! 
    544                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2  
     544                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
    545545                     &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    546546                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
     
    565565               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   & 
    566566               &               / fse3w(:,:,jk) * tmask(:,:,jk) 
    567          END DO  
     567         END DO 
    568568#if defined key_zdfddm 
    569569         DO jk = 2, jpkm1                                 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 
    570570            DO jj = 1, jpj 
    571571               DO ji = 1, jpi 
    572                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )   
     572                  zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 
    573573                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 
    574574                  rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds 
     
    587587 
    588588 
    589    SUBROUTINE eos_alpbet( pts, palph, pbeta ) 
    590       !!---------------------------------------------------------------------- 
    591       !!                 ***  ROUTINE ldf_slp_grif  *** 
    592       !! 
    593       !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points.  
    594       !! 
    595       !! ** Method  :   calculates alpha and beta at T-points  
     589   SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 
     590      !!---------------------------------------------------------------------- 
     591      !!                 ***  ROUTINE eos_alpbet  *** 
     592      !! 
     593      !! ** Purpose :   Calculates the in situ thermal/haline expansion ratio at T-points 
     594      !! 
     595      !! ** Method  :   calculates alpha / beta ratio at T-points 
    596596      !!       * nn_eos = 0  : UNESCO sea water properties 
    597       !!         The brunt-vaisala frequency is computed using the polynomial 
    598       !!      polynomial expression of McDougall (1987): 
    599       !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 
    600       !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 
    601       !!      computed and used in zdfddm module : 
    602       !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 
     597      !!                       The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 
     598      !!                       polynomial expression of McDougall (1987). 
     599      !!                       Scalar beta0 is returned = 1. 
    603600      !!       * nn_eos = 1  : linear equation of state (temperature only) 
    604       !!            N^2 = grav * rn_alpha * dk[ t ]/e3w 
     601      !!                       The ratio is undefined, so we return alpha as palpbet 
     602      !!                       Scalar beta0 is returned = 0. 
    605603      !!       * nn_eos = 2  : linear equation of state (temperature & salinity) 
    606       !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 
    607       !!       * nn_eos = 3  : Jackett JAOT 2003 ??? 
    608       !! 
    609       !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point 
    610       !!---------------------------------------------------------------------- 
    611       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity 
    612       REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff. 
    613       ! 
     604      !!                       The alpha/beta ratio is returned as ralpbet 
     605      !!                       Scalar beta0 is returned = 1. 
     606      !! 
     607      !! ** Action  : - palpbet : thermal/haline expansion ratio at T-points 
     608      !!            :   beta0   : 1. or 0. 
     609      !!---------------------------------------------------------------------- 
     610      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts       ! pot. temperature & salinity 
     611      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palpbet   ! thermal/haline expansion ratio 
     612      REAL(wp),                              INTENT(  out) ::   beta0     ! set = 1 except with case 1 eos, rho=rho(T) 
     613      !! 
    614614      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    615       REAL(wp) ::   zt, zs, zh   ! local scalars  
     615      REAL(wp) ::   zt, zs, zh   ! local scalars 
    616616      !!---------------------------------------------------------------------- 
    617617      ! 
     
    624624                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature 
    625625                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35) 
    626                   zh = fsdept(ji,jj,jk)              ! depth in meters  
    627                   ! 
    628                   pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   & 
    629                      &                                        - 0.301985e-05_wp ) * zt   & 
    630                      &           + 0.785567e-03_wp                                       & 
    631                      &           + (     0.515032e-08_wp * zs                            & 
    632                      &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   & 
    633                      &           + ( (   0.121551e-17_wp * zh                            & 
    634                      &                 - 0.602281e-15_wp * zs                            & 
    635                      &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   & 
    636                      &                                        + 0.408195e-10_wp   * zs   & 
    637                      &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   & 
    638                      &                                        - 0.121555e-07_wp ) * zh 
    639                      ! 
    640                   palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             & 
    641                       &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
    642                       &                                  - 0.203814e-03_wp ) * zt   & 
    643                       &                                  + 0.170907e-01_wp ) * zt   & 
    644                       &   + 0.665157e-01_wp                                         & 
    645                       &   +     ( - 0.678662e-05_wp * zs                            & 
    646                       &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
    647                       &   +   ( ( - 0.302285e-13_wp * zh                            & 
    648                       &           - 0.251520e-11_wp * zs                            & 
    649                       &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
    650                       &           - 0.164759e-06_wp * zs                            & 
    651                       &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
    652                       &                                  + 0.380374e-04_wp ) * zh) 
     626                  zh = fsdept(ji,jj,jk)               ! depth in meters 
     627                  ! 
     628                  palpbet(ji,jj,jk) =                                              & 
     629                     &     ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   & 
     630                     &                                  - 0.203814e-03_wp ) * zt   & 
     631                     &                                  + 0.170907e-01_wp ) * zt   & 
     632                     &   + 0.665157e-01_wp                                         & 
     633                     &   +     ( - 0.678662e-05_wp * zs                            & 
     634                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   & 
     635                     &   +   ( ( - 0.302285e-13_wp * zh                            & 
     636                     &           - 0.251520e-11_wp * zs                            & 
     637                     &           + 0.512857e-12_wp * zt * zt              ) * zh   & 
     638                     &           - 0.164759e-06_wp * zs                            & 
     639                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   & 
     640                     &                                  + 0.380374e-04_wp ) * zh 
    653641               END DO 
    654642            END DO 
    655643         END DO 
    656          ! 
    657       CASE ( 1 ) 
    658          palph(:,:,:) = - rn_alpha 
    659          pbeta(:,:,:) =   0._wp 
    660          ! 
    661       CASE ( 2 ) 
    662          palph(:,:,:) = - rn_alpha 
    663          pbeta(:,:,:) =   rn_beta 
     644         beta0 = 1._wp 
     645         ! 
     646      CASE ( 1 )              !==  Linear formulation = F( temperature )  ==! 
     647         palpbet(:,:,:) = rn_alpha 
     648         beta0 = 0._wp 
     649         ! 
     650      CASE ( 2 )              !==  Linear formulation = F( temperature , salinity )  ==! 
     651         palpbet(:,:,:) = ralpbet 
     652         beta0 = 1._wp 
    664653         ! 
    665654      CASE DEFAULT 
     
    747736 
    748737   !!====================================================================== 
    749 END MODULE eosbn2   
     738END MODULE eosbn2 
  • branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90

    r2715 r2840  
    44   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend 
    55   !!====================================================================== 
    6    !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)   
     6   !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec) 
    77   !!                !          Griffies operator version adapted from traldf_iso.F90 
    88   !!---------------------------------------------------------------------- 
     
    1111   !!   'key_ldfslp'               slope of the lateral diffusive direction 
    1212   !!---------------------------------------------------------------------- 
    13    !!   tra_ldf_iso_grif  : update the tracer trend with the horizontal component   
    14    !!                       of the Griffies iso-neutral laplacian operator  
     13   !!   tra_ldf_iso_grif  : update the tracer trend with the horizontal component 
     14   !!                       of the Griffies iso-neutral laplacian operator 
    1515   !!---------------------------------------------------------------------- 
    1616   USE oce             ! ocean dynamics and active tracers 
     
    5353      !!                  ***  ROUTINE tra_ldf_iso_grif  *** 
    5454      !! 
    55       !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
    56       !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and  
     55      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
     56      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 
    5757      !!      add it to the general trend of tracer equation. 
    5858      !! 
    59       !! ** Method  :   The horizontal component of the lateral diffusive trends  
     59      !! ** Method  :   The horizontal component of the lateral diffusive trends 
    6060      !!      is provided by a 2nd order operator rotated along neural or geopo- 
    6161      !!      tential surfaces to which an eddy induced advection can be added 
     
    6767      !! 
    6868      !!      2nd part :  horizontal fluxes of the lateral mixing operator 
    69       !!      ========     
     69      !!      ======== 
    7070      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 
    7171      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ] 
     
    9999      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    100100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    101       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend 
    102102      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
    103103      ! 
     
    143143 
    144144      !!---------------------------------------------------------------------- 
    145       !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv  
     145      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv 
    146146      !!---------------------------------------------------------------------- 
    147147 
     
    159159               DO jj = 1, jpjm1 
    160160                  DO ji = 1, fs_jpim1 
     161                     ze1ur = 1._wp / e1u(ji,jj) 
    161162                     ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
    162163                     zbu   = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    163164                     zah   = fsahtu(ji,jj,jk)                                  !  fsaht(ji+ip,jj,jk) 
    164165                     zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    165                      zslope2 = zslope_skew - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 
     166                     ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
     167                     !    (do this by *adding* gradient of depth) 
     168                     zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 
    166169                     zslope2 = zslope2 *zslope2 
    167170                     ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp)    & 
     
    182185               DO jj = 1, jpjm1 
    183186                  DO ji=1,fs_jpim1 
     187                     ze2vr = 1._wp / e2v(ji,jj) 
    184188                     ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 
    185189                     zbv   = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    186                      zah   = fsahtu(ji,jj,jk)                                       !fsaht(ji,jj+jp,jk) 
     190                     zah   = fsahtv(ji,jj,jk)                                       !fsaht(ji,jj+jp,jk) 
    187191                     zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    188                      zslope2 = zslope_skew - ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 
     192                     ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
     193                     !    (do this by *adding* gradient of depth) 
     194                     zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 
    189195                     zslope2 = zslope2 * zslope2 
    190196                     ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp)   & 
     
    207213         zftu(:,:,:) = 0._wp 
    208214         zftv(:,:,:) = 0._wp 
    209          !                                                
     215         ! 
    210216         DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==! 
    211217            DO jj = 1, jpjm1 
     
    216222            END DO 
    217223         END DO 
    218          IF( ln_zps ) THEN                               ! partial steps: correction at the last level 
     224         IF( ln_zps.and.ln_grad_zps ) THEN              ! partial steps: correction at the last level 
    219225# if defined key_vectopt_loop 
    220226            DO jj = 1, 1 
     
    224230               DO ji = 1, jpim1 
    225231# endif 
    226                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    227                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)       
     232                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
     233                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    228234               END DO 
    229235            END DO 
     
    244250            ENDIF 
    245251 
    246             IF( .NOT. l_triad_iso ) THEN 
    247                triadi = triadi_g 
    248                triadj = triadj_g 
    249             ENDIF 
    250  
    251             DO ip = 0, 1              !==  Horizontal & vertical fluxes 
    252                DO kp = 0, 1 
    253                   DO jj = 1, jpjm1 
    254                      DO ji = 1, fs_jpim1 
    255                         ze1ur = 1._wp / e1u(ji,jj) 
    256                         zdxt = zdit(ji,jj,jk) * ze1ur 
    257                         ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
    258                         zdzt  = zdkt(ji+ip,jj,kp) * ze3wr 
    259                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    260                         zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
    261  
    262                         zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    263                         zah = fsahtu(ji,jj,jk)   !*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ???? 
    264                         zah_slp  = zah * zslope_iso 
    265                         zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(ji+ip,jj,jk)*zslope_skew 
    266                         zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    267                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
     252 
     253            IF (ln_botmix) THEN 
     254               DO ip = 0, 1              !==  Horizontal & vertical fluxes 
     255                  DO kp = 0, 1 
     256                     DO jj = 1, jpjm1 
     257                        DO ji = 1, fs_jpim1 
     258                           ze1ur = 1._wp / e1u(ji,jj) 
     259                           zdxt = zdit(ji,jj,jk) * ze1ur 
     260                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
     261                           zdzt  = zdkt(ji+ip,jj,kp) * ze3wr 
     262                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     263                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
     264 
     265                           zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     266                           ! ln_botmix is .T. don't mask zah for bottom half cells 
     267                           zah = fsahtu(ji,jj,jk)   !*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ???? 
     268                           zah_slp  = zah * zslope_iso 
     269                           zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(ji+ip,jj,jk)*zslope_skew 
     270                           zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
     271                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
     272                        END DO 
    268273                     END DO 
    269274                  END DO 
    270275               END DO 
    271             END DO 
    272  
    273             DO jp = 0, 1 
    274                DO kp = 0, 1 
    275                   DO jj = 1, jpjm1 
    276                      DO ji = 1, fs_jpim1 
    277                         ze2vr = 1._wp / e2v(ji,jj) 
    278                         zdyt = zdjt(ji,jj,jk) * ze2vr 
    279                         ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 
    280                         zdzt = zdkt(ji,jj+jp,kp) * ze3wr 
    281                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    282                         zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    283                         zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    284                         zah = fsahtv(ji,jj,jk)        !*vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk) 
    285                         zah_slp = zah * zslope_iso 
    286                         zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew    !fsaeit(ji,jj+jp,jk)*zslope_skew 
    287                         zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    288                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
     276 
     277               DO jp = 0, 1 
     278                  DO kp = 0, 1 
     279                     DO jj = 1, jpjm1 
     280                        DO ji = 1, fs_jpim1 
     281                           ze2vr = 1._wp / e2v(ji,jj) 
     282                           zdyt = zdjt(ji,jj,jk) * ze2vr 
     283                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 
     284                           zdzt = zdkt(ji,jj+jp,kp) * ze3wr 
     285                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     286                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
     287                           zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     288                           ! ln_botmix is .T. don't mask zah for bottom half cells 
     289                           zah = fsahtv(ji,jj,jk)        !*vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk) 
     290                           zah_slp = zah * zslope_iso 
     291                           zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew    !fsaeit(ji,jj+jp,jk)*zslope_skew 
     292                           zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
     293                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
     294                        END DO 
    289295                     END DO 
    290296                  END DO 
    291297               END DO 
    292             END DO 
    293  
     298            ELSE 
     299               DO ip = 0, 1              !==  Horizontal & vertical fluxes 
     300                  DO kp = 0, 1 
     301                     DO jj = 1, jpjm1 
     302                        DO ji = 1, fs_jpim1 
     303                           ze1ur = 1._wp / e1u(ji,jj) 
     304                           zdxt = zdit(ji,jj,jk) * ze1ur 
     305                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 
     306                           zdzt  = zdkt(ji+ip,jj,kp) * ze3wr 
     307                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     308                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
     309 
     310                           zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     311                           ! ln_botmix is .F. mask zah for bottom half cells 
     312                           zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ???? 
     313                           zah_slp  = zah * zslope_iso 
     314                           zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(ji+ip,jj,jk)*zslope_skew 
     315                           zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
     316                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
     317                        END DO 
     318                     END DO 
     319                  END DO 
     320               END DO 
     321 
     322               DO jp = 0, 1 
     323                  DO kp = 0, 1 
     324                     DO jj = 1, jpjm1 
     325                        DO ji = 1, fs_jpim1 
     326                           ze2vr = 1._wp / e2v(ji,jj) 
     327                           zdyt = zdjt(ji,jj,jk) * ze2vr 
     328                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 
     329                           zdzt = zdkt(ji,jj+jp,kp) * ze3wr 
     330                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     331                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
     332                           zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     333                           ! ln_botmix is .F. mask zah for bottom half cells 
     334                           zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk) 
     335                           zah_slp = zah * zslope_iso 
     336                           zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew    !fsaeit(ji,jj+jp,jk)*zslope_skew 
     337                           zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
     338                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
     339                        END DO 
     340                     END DO 
     341                  END DO 
     342               END DO 
     343            END IF 
    294344            !                        !==  divergence and add to the general trend  ==! 
    295345            DO jj = 2 , jpjm1 
     
    320370#if defined key_diaar5 
    321371         IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    322             z2d(:,:) = 0._wp  
    323             zztmp = rau0 * rcp  
     372            z2d(:,:) = 0._wp 
     373            zztmp = rau0 * rcp 
    324374            DO jk = 1, jpkm1 
    325375               DO jj = 2, jpjm1 
    326376                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    327                      z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     377                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
    328378                  END DO 
    329379               END DO 
     
    332382            CALL lbc_lnk( z2d, 'U', -1. ) 
    333383            CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    334             z2d(:,:) = 0._wp  
     384            z2d(:,:) = 0._wp 
    335385            DO jk = 1, jpkm1 
    336386               DO jj = 2, jpjm1 
    337387                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    338                      z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     388                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
    339389                  END DO 
    340390               END DO 
  • branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90

    r2715 r2840  
    7676         WRITE(charout, FMT="('ldf0 ')") ;  CALL prt_ctl_trc_info(charout) 
    7777                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 
    78          CALL tra_ldf_iso   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 
     78         IF( ln_traldf_grif ) THEN 
     79            CALL tra_ldf_iso_grif( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 
     80         ELSE 
     81            CALL tra_ldf_iso   ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 
     82         ENDIF 
    7983         WRITE(charout, FMT="('ldf1 ')") ;  CALL prt_ctl_trc_info(charout) 
    8084                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 
Note: See TracChangeset for help on using the changeset viewer.