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 2450 for branches/nemo_v3_3_beta – NEMO

Ignore:
Timestamp:
2010-12-04T16:20:50+01:00 (13 years ago)
Author:
gm
Message:

v3.3beta: #766 share the deepest ocean level indices continuaton

Location:
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC
Files:
23 edited

Legend:

Unmodified
Added
Removed

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r2287 r2450  
    1616   USE oce            ! ocean dynamics and active tracers  
    1717   USE dom_oce        ! ocean space and time domain 
    18    USE eosbn2          ! equation of state                (eos_bn2 routine) 
     18   USE eosbn2         ! equation of state                (eos_bn2 routine) 
    1919   USE lib_mpp        ! distribued memory computing library 
    2020   USE iom            ! I/O manager library 
    … …  
    2828   LOGICAL, PUBLIC, PARAMETER :: lk_diaar5 = .TRUE.   ! coupled flag 
    2929 
    30    REAL(wp)                         ::   vol0               ! ocean volume (interior domain) 
    31    REAL(wp)                         ::   area_tot           ! total ocean surface (interior domain) 
    32    REAL(wp), DIMENSION(jpi,jpj    ) ::   area               ! cell surface (interior domain) 
    33    REAL(wp), DIMENSION(jpi,jpj    ) ::   thick0             ! ocean thickness (interior domain) 
    34    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sn0                ! initial salinity 
     30   REAL(wp)                         ::   vol0         ! ocean volume (interior domain) 
     31   REAL(wp)                         ::   area_tot     ! total ocean surface (interior domain) 
     32   REAL(wp), DIMENSION(jpi,jpj    ) ::   area         ! cell surface (interior domain) 
     33   REAL(wp), DIMENSION(jpi,jpj    ) ::   thick0       ! ocean thickness (interior domain) 
     34   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sn0          ! initial salinity 
    3535       
    3636   !! * Substitutions 
    … …  
    3939   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4040   !! $Id$ 
    41    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    42    !!---------------------------------------------------------------------- 
    43  
     41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     42   !!---------------------------------------------------------------------- 
    4443CONTAINS 
    4544 
    … …  
    7776      CALL eos( ztsn, zrhd )                       ! now in situ density using initial salinity 
    7877      ! 
    79       zbotpres(:,:) = 0.e0                            ! no atmospheric surface pressure, levitating sea-ice 
     78      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    8079      DO jk = 1, jpkm1 
    8180         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
    8281      END DO 
    83       IF( .NOT. lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     82      IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
    8483      !                                          
    8584      zarho = SUM( area(:,:) * zbotpres(:,:) )  
    … …  
    9089      !                                         ! steric sea surface height 
    9190      CALL eos( tsn, zrhd, zrhop )                 ! now in situ and potential density 
    92       zrhop(:,:,jpk) = 0.e0 
     91      zrhop(:,:,jpk) = 0._wp 
    9392      CALL iom_put( 'rhop', zrhop ) 
    9493      ! 
    95       zbotpres(:,:) = 0.e0                            ! no atmospheric surface pressure, levitating sea-ice 
     94      zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice 
    9695      DO jk = 1, jpkm1 
    9796         zbotpres(:,:) = zbotpres(:,:) + fse3t(:,:,jk) * zrhd(:,:,jk) 
    9897      END DO 
    99       IF( .NOT. lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     98      IF( .NOT.lk_vvl )   zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
    10099      !     
    101100      zarho = SUM( area(:,:) * zbotpres(:,:) )  
    … …  
    105104       
    106105      !                                         ! ocean bottom pressure 
    107       zztmp = rau0 * grav * 1.e-4                     ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
     106      zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
    108107      zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 
    109108      CALL iom_put( 'botpres', zbotpres ) 
    110109 
    111110      !                                         ! Mean density anomalie, temperature and salinity 
    112       ztemp = 0.e0 
    113       zsal  = 0.e0 
     111      ztemp = 0._wp 
     112      zsal  = 0._wp 
    114113      DO jk = 1, jpkm1 
    115114         DO jj = 1, jpj 
    116115            DO ji = 1, jpi 
    117116               zztmp = area(ji,jj) * fse3t(ji,jj,jk) 
    118                ztemp = ztemp + zztmp * tn  (ji,jj,jk) 
    119                zsal  = zsal  + zztmp * sn  (ji,jj,jk) 
     117               ztemp = ztemp + zztmp * tn(ji,jj,jk) 
     118               zsal  = zsal  + zztmp * sn(ji,jj,jk) 
    120119            END DO 
    121120         END DO 
    122121      END DO 
    123       IF( .NOT. lk_vvl ) THEN 
    124          ztemp = ztemp + SUM( zarea_ssh(:,:) * tn  (:,:,1) ) 
    125          zsal  = zsal  + SUM( zarea_ssh(:,:) * sn  (:,:,1) ) 
     122      IF( .NOT.lk_vvl ) THEN 
     123         ztemp = ztemp + SUM( zarea_ssh(:,:) * tn(:,:,1) ) 
     124         zsal  = zsal  + SUM( zarea_ssh(:,:) * sn(:,:,1) ) 
    126125      ENDIF 
    127126      IF( lk_mpp ) THEN   
    … …  
    146145      !!                    
    147146      !! ** Purpose :   initialization for AR5 diagnostic computation 
    148       !! 
    149147      !!---------------------------------------------------------------------- 
    150148      INTEGER  ::   inum 
    … …  
    159157      area_tot = SUM( area(:,:) )   ;   IF( lk_mpp )   CALL mpp_sum( area_tot ) 
    160158 
    161       vol0        = 0.e0 
    162       thick0(:,:) = 0.e0 
     159      vol0        = 0._wp 
     160      thick0(:,:) = 0._wp 
    163161      DO jk = 1, jpkm1 
    164162         vol0        = vol0        + SUM( area (:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk) ) 
    … …  
    171169      CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,2), 12 ) 
    172170      CALL iom_close( inum ) 
    173       sn0(:,:,:) = 0.5 * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
    174       sn0(:,:,:) = sn0(:,:,:)*tmask(:,:,:) 
     171      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
     172      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:) 
    175173      IF( ln_zps ) THEN               ! z-coord. partial steps 
    176174         DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    177175            DO ji = 1, jpi 
    178                ik = mbathy(ji,jj) - 1 
    179                IF( ik > 2 ) THEN 
     176               ik = mbkt(ji,jj) 
     177               IF( ik > 1 ) THEN 
    180178                  zztmp = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
    181                   sn0(ji,jj,ik) = (1.-zztmp) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
     179                  sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
    182180               ENDIF 
    183181            END DO 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

    r2287 r2450  
    1111   !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag 
    1212   !!---------------------------------------------------------------------- 
    13  
    1413#if   defined key_diahth   ||   defined key_esopa 
    1514   !!---------------------------------------------------------------------- 
    … …  
    1817   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer 
    1918   !!---------------------------------------------------------------------- 
    20    !! * Modules used 
    2119   USE oce             ! ocean dynamics and tracers 
    2220   USE dom_oce         ! ocean space and time domain 
    2321   USE phycst          ! physical constants 
    2422   USE in_out_manager  ! I/O manager 
    25    USE iom 
     23   USE iom             ! I/O library 
    2624 
    2725   IMPLICIT NONE 
    2826   PRIVATE 
    2927 
    30    !! * Routine accessibility 
    31    PUBLIC dia_hth    ! routine called by step.F90 
    32  
    33    !! * Shared module variables 
     28   PUBLIC   dia_hth    ! routine called by step.F90 
     29 
    3430   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag 
    3531   ! note: following variables should move to local variables once iom_put is always used  
    … …  
    4440   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    4541   !! $Id$  
    46    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    47    !!---------------------------------------------------------------------- 
    48  
     42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     43   !!---------------------------------------------------------------------- 
    4944CONTAINS 
    5045 
    … …  
    6863      !! 
    6964      !! ** Method :  
    70       !! 
    7165      !!------------------------------------------------------------------- 
    7266      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    7367      !! 
    7468      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments 
    75       INTEGER                          ::   iid, iif, ilevel      ! temporary integers 
     69      INTEGER                          ::   iid, ilevel           ! temporary integers 
    7670      INTEGER, DIMENSION(jpi,jpj)      ::   ik20, ik28            ! levels 
    7771      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth 
    … …  
    10397 
    10498      ! initialization 
    105       ztinv  (:,:) = 0.e0_wp   
    106       zdepinv(:,:) = 0.e0_wp   
    107       zmaxdzT(:,:) = 0.e0_wp   
     99      ztinv  (:,:) = 0._wp   
     100      zdepinv(:,:) = 0._wp   
     101      zmaxdzT(:,:) = 0._wp   
    108102      DO jj = 1, jpj 
    109103         DO ji = 1, jpi 
    … …  
    128122      ! Preliminary computation 
    129123      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC) 
    130       DO jj=1, jpj 
    131          DO ji=1, jpi 
     124      DO jj = 1, jpj 
     125         DO ji = 1, jpi 
    132126            IF( tmask(ji,jj,nla10) == 1. ) THEN 
    133127               zu  =  1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10)   & 
    … …  
    139133               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw) 
    140134            ELSE 
    141                zdelr(ji,jj) = 0.e0 
     135               zdelr(ji,jj) = 0._wp 
    142136            ENDIF 
    143137         END DO 
    … …  
    153147         DO jj = 1, jpj 
    154148            DO ji = 1, jpi 
    155  
     149               ! 
    156150               zzdep = fsdepw(ji,jj,jk) 
    157151               zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz) 
    … …  
    189183         DO jj = 1, jpj 
    190184            DO ji = 1, jpi 
    191  
     185               ! 
    192186               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1) 
    193  
     187               ! 
    194188               zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk)                  ! - delta T(10m) 
    195189               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2 
    … …  
    203197               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03 
    204198               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2 
    205  
     199               ! 
    206200            END DO 
    207201         END DO 
    … …  
    237231      DO jj = 1, jpj 
    238232         DO ji = 1, jpi 
    239  
    240             iif = mbathy(ji,jj) 
    241             zzdep = fsdepw(ji,jj,iif) 
    242  
     233            ! 
     234            zzdep = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom 
     235            ! 
    243236            iid = ik20(ji,jj) 
    244237            IF( iid /= 1 ) THEN  
    245                ! linear interpolation 
    246                zztmp =      fsdept(ji,jj,iid  )   & 
     238               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation 
    247239                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   & 
    248240                  &  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   & 
    249241                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 
    250                ! bound by the ocean depth, minimum value, first T-point depth 
    251                hd20(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep) 
     242               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth 
    252243            ELSE  
    253                hd20(ji,jj)=0. 
     244               hd20(ji,jj) = 0._wp 
    254245            ENDIF 
    255  
     246            ! 
    256247            iid = ik28(ji,jj) 
    257248            IF( iid /= 1 ) THEN  
    258                ! linear interpolation 
    259                zztmp =      fsdept(ji,jj,iid  )   & 
     249               zztmp =      fsdept(ji,jj,iid  )   &                     ! linear interpolation 
    260250                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   & 
    261251                  &  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   & 
    262252                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) ) 
    263                ! bound by the ocean depth, minimum value, first T-point depth 
    264                hd28(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep ) 
     253               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth 
    265254            ELSE  
    266                hd28(ji,jj) = 0. 
     255               hd28(ji,jj) = 0._wp 
    267256            ENDIF 
    268257 
    … …  
    277266 
    278267      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...) 
    279       ilevel = 0 
    280       zthick_0 = 0.e0_wp 
     268      ilevel   = 0 
     269      zthick_0 = 0._wp 
    281270      DO jk = 1, jpkm1                       
    282271         zthick_0 = zthick_0 + e3t_0(jk) 
    … …  
    284273      END DO 
    285274      ! surface boundary condition 
    286       IF( lk_vvl ) THEN   ;   zthick(:,:) = 0.e0_wp     ;   htc3(:,:) = 0.e0_wp                                    
     275      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                    
    287276      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)    
    288277      ENDIF 
    … …  
    303292      htc3(:,:) = zcoef * htc3(:,:) 
    304293      CALL iom_put( "hc300", htc3 )      ! first 300m heat content 
    305  
    306  
     294      ! 
    307295   END SUBROUTINE dia_hth 
    308296 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r2287 r2450  
    44   !! Ocean initialization : write the ocean domain mesh ask file(s) 
    55   !!====================================================================== 
     6   !! History :  OPA  ! 1997-02  (G. Madec)  Original code 
     7   !!            8.1  ! 1999-11  (M. Imbard)  NetCDF FORMAT with IOIPSL 
     8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file 
     9   !!---------------------------------------------------------------------- 
    610 
    711   !!---------------------------------------------------------------------- 
    … …  
    1115   !!                         = 3  :   mesh_hgr, mesh_zgr and mask 
    1216   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1417   USE dom_oce         ! ocean space and time domain 
    15    USE in_out_manager 
    16    USE iom 
    17    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    18    USE lib_mpp 
     18   USE in_out_manager  ! I/O manager 
     19   USE iom             ! I/O library 
     20   USE lbclnk          ! lateral boundary conditions - mpp exchanges 
     21   USE lib_mpp         ! MPP library 
    1922 
    2023   IMPLICIT NONE 
    … …  
    2831   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    2932   !! $Id$  
    30    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    31    !!---------------------------------------------------------------------- 
    32  
     33   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     34   !!---------------------------------------------------------------------- 
    3335CONTAINS 
    3436 
    … …  
    5860      !!                        thickness of the bottom points hdep[tw] and e3[tw]_ps 
    5961      !! 
    60       !! ** output file :  
    61       !!      meshmask.nc  : domain size, horizontal grid-point position, 
    62       !!                     masks, depth and vertical scale factors 
    63       !! 
    64       !! History : 
    65       !!        !  97-02  (G. Madec)  Original code 
    66       !!        !  99-11  (M. Imbard)  NetCDF FORMAT with IOIPSL 
    67       !!   9.0  !  02-08  (G. Madec)  F90 and several file 
     62      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position, 
     63      !!                                   masks, depth and vertical scale factors 
    6864      !!---------------------------------------------------------------------- 
    6965      INTEGER                          ::   inum0    ! temprary units for 'mesh_mask.nc' file 
    … …  
    7268      INTEGER                          ::   inum3    ! temprary units for 'mesh_hgr.nc'  file 
    7369      INTEGER                          ::   inum4    ! temprary units for 'mesh_zgr.nc'  file 
    74       INTEGER                          ::   ji, jj, jk, ik 
     70      INTEGER                          ::   ji, jj, jk   ! dummy loop indices 
    7571      REAL(wp), DIMENSION(jpi,jpj)     ::   zprt     ! temporary array for bathymetry  
    7672      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu    ! 3D depth of U point 
    … …  
    118114         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib ) 
    119115         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib ) 
    120           
     116         ! 
    121117      END SELECT 
    122118       
    … …  
    160156      CALL iom_rstput( 0, 0, inum3, 'ff', ff, ktype = jp_r8 )           !    ! coriolis factor 
    161157       
    162 !       note that mbathy has been modified in dommsk or in solver. 
    163 !       it is the number of non-zero "w" levels in the water, and the minimum  
    164 !       value (on land) is 2. We define zprt as the number of "T" points in the ocean  
    165 !       at any location, and zero on land.  
    166 ! 
    167       zprt = tmask(:,:,1)*(mbathy-1) 
    168       CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 ) 
     158      ! note that mbkt is set to 1 over land ==> use surface tmask 
     159      zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp ) 
     160      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points 
    169161             
    170162      IF( ln_sco ) THEN                                         ! s-coordinate 
    … …  
    173165         CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 
    174166         CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) 
    175           
     167         ! 
    176168         CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt )         !    ! scaling coef. 
    177169         CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )   
    … …  
    179171         CALL iom_rstput( 0, 0, inum4, 'esigt', esigt ) 
    180172         CALL iom_rstput( 0, 0, inum4, 'esigw', esigw ) 
    181           
     173         ! 
    182174         CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )             !    ! scale factors 
    183175         CALL iom_rstput( 0, 0, inum4, 'e3u', e3u ) 
    184176         CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 
    185177         CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    186           
     178         ! 
    187179         CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 )    !    ! stretched system 
    188180         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 
    … …  
    190182       
    191183      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps 
    192  
     184         ! 
    193185         IF( nmsh <= 6 ) THEN                                   !    ! 3D vertical scale factors 
    194186            CALL iom_rstput( 0, 0, inum4, 'e3t', e3t )          
    … …  
    196188            CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 
    197189            CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 
    198          ELSE                                                   !    ! 2D bottom scale factors 
    199             DO jj = 1,jpj   ;   DO ji = 1,jpi 
    200                ik = NINT( zprt(ji,jj) )   ! take care that mbathy is not what you think it is here ! 
    201                IF ( ik /= 0 ) THEN   ;   e3tp(ji,jj) = e3t(ji,jj,ik)   ;   e3wp(ji,jj) = e3w(ji,jj,ik) 
    202                ELSE                  ;   e3tp(ji,jj) = 0.              ;   e3wp(ji,jj) = 0. 
    203                ENDIF 
    204             END DO   ;   END DO 
     190         ELSE                                                   !    ! 2D masked bottom ocean scale factors 
     191            DO jj = 1,jpj    
     192               DO ji = 1,jpi 
     193                  e3tp(ji,jj) = e3t(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
     194                  e3wp(ji,jj) = e3w(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
     195               END DO 
     196            END DO 
    205197            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )       
    206198            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp ) 
    207199         END IF 
    208  
     200         ! 
    209201         IF( nmsh <= 3 ) THEN                                   !    ! 3D depth 
    210202            CALL iom_rstput( 0, 0, inum4, 'gdept', gdept, ktype = jp_r4 )      
    211             DO jk = 1,jpk   ;   DO jj = 1, jpjm1   ;   DO ji = 1, fs_jpim1   ! vector opt. 
    212                zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji+1,jj  ,jk) ) 
    213                zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk), gdept(ji  ,jj+1,jk) ) 
    214             END DO   ;   END DO   ;   END DO 
     203            DO jk = 1,jpk    
     204               DO jj = 1, jpjm1    
     205                  DO ji = 1, fs_jpim1   ! vector opt. 
     206                     zdepu(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji+1,jj  ,jk) ) 
     207                     zdepv(ji,jj,jk) = MIN( gdept(ji,jj,jk) , gdept(ji  ,jj+1,jk) ) 
     208                  END DO    
     209               END DO    
     210            END DO 
    215211            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. )  
    216212            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r4 ) 
    … …  
    218214            CALL iom_rstput( 0, 0, inum4, 'gdepw', gdepw, ktype = jp_r4 ) 
    219215         ELSE                                                   !    ! 2D bottom depth 
    220             DO jj = 1,jpj   ;   DO ji = 1,jpi 
    221                ik = NINT( zprt(ji,jj) )   ! take care that mbathy is not what you think it is here ! 
    222                IF ( ik /= 0 ) THEN   ;   hdept(ji,jj) = gdept(ji,jj,ik)   ;   hdepw(ji,jj) = gdepw(ji,jj,ik+1) 
    223                ELSE                  ;   hdept(ji,jj) = 0.                ;   hdepw(ji,jj) = 0. 
    224                ENDIF 
    225             END DO   ;   END DO 
     216            DO jj = 1,jpj    
     217               DO ji = 1,jpi 
     218                  hdept(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1) 
     219                  hdepw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1) 
     220               END DO 
     221            END DO 
    226222            CALL iom_rstput( 0, 0, inum4, 'hdept' , hdept, ktype = jp_r4 )      
    227223            CALL iom_rstput( 0, 0, inum4, 'hdepw' , hdepw, ktype = jp_r4 )  
    228224         ENDIF 
    229  
     225         ! 
    230226         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0 )     !    ! reference z-coord. 
    231227         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0 ) 
    … …  
    233229         CALL iom_rstput( 0, 0, inum4, 'e3w_0'  , e3w_0   ) 
    234230      ENDIF 
    235        
    236231       
    237232      IF( ln_zco ) THEN 
    … …  
    256251         CALL iom_close( inum4 ) 
    257252      END SELECT 
    258        
     253      ! 
    259254   END SUBROUTINE dom_wri 
    260255 
    … …  
    268263      !! ** Method  :   1) aplly lbc_lnk on an array with different values for each element 
    269264      !!                2) check which elements have been changed 
    270       !! 
    271265      !!---------------------------------------------------------------------- 
    272266      CHARACTER(len=1)            , INTENT(in   ) ::  cdgrd   !  
    273267      REAL(wp), DIMENSION(jpi,jpj)                ::  puniq   !  
    274  
     268      ! 
    275269      REAL(wp), DIMENSION(jpi,jpj  ) ::  ztstref   ! array with different values for each element  
    276270      REAL(wp)                       ::  zshift    ! shift value link to the process number 
    … …  
    278272      INTEGER                        ::  ji        ! dummy loop indices 
    279273      !!---------------------------------------------------------------------- 
    280  
     274      ! 
    281275      ! build an array with different values for each element  
    282276      ! in mpp: make sure that these values are different even between process 
    283277      ! -> apply a shift value according to the process number 
    284278      zshift = jpi * jpj * ( narea - 1 ) 
    285       ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji, wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) ) 
    286     
     279      ztstref(:,:) = RESHAPE( (/ (zshift + REAL(ji,wp), ji = 1, jpi*jpj) /), (/ jpi, jpj /) ) 
     280      ! 
    287281      puniq(:,:) = ztstref(:,:)                   ! default definition 
    288282      CALL lbc_lnk( puniq, cdgrd, 1. )            ! apply boundary conditions 
    289283      lldbl(:,:,1) = puniq(:,:) == ztstref(:,:)   ! check which values have been changed  
    290        
     284      ! 
    291285      puniq(:,:) = 1.                             ! default definition 
    292286      ! fill only the inner part of the cpu with llbl converted into real  
    293       puniq(nldi:nlei,nldj:nlej) = REAL(COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ), wp) 
    294  
     287      puniq(nldi:nlei,nldj:nlej) = REAL( COUNT( lldbl(nldi:nlei,nldj:nlej,:), dim = 3 ) , wp ) 
     288      ! 
    295289   END FUNCTION dom_uniq 
    296  
    297290 
    298291   !!====================================================================== 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtasal.F90

    r2436 r2450  
    190190            DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step) 
    191191               DO ji = 1, jpi 
    192                   ik = mbathy(ji,jj) - 1 
    193                   IF( ik > 2 ) THEN 
     192                  ik = mbkt(ji,jj) 
     193                  IF( ik > 1 ) THEN 
    194194                     zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
    195195                     s_dta(ji,jj,ik) = (1.-zl) * s_dta(ji,jj,ik) + zl * s_dta(ji,jj,ik-1) 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtatem.F90

    r2436 r2450  
    203203            DO jj = 1, jpj                ! interpolation of temperature at the last level 
    204204               DO ji = 1, jpi 
    205                   ik = mbathy(ji,jj) - 1 
    206                   IF( ik > 2 ) THEN 
     205                  ik = mbkt(ji,jj) 
     206                  IF( ik > 1 ) THEN 
    207207                     zl = ( gdept_0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) ) 
    208208                     t_dta(ji,jj,ik) = (1.-zl) * t_dta(ji,jj,ik) + zl * t_dta(ji,jj,ik-1) 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r2287 r2450  
    1111   !!            8.5  !  2002-08  (A. Bozec)  hpg_zps: Original code 
    1212   !!   NEMO     1.0  !  2005-10  (A. Beckmann, B.W. An)  various s-coordinate options 
    13    !!                           Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot  
     13   !!                 !         Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot  
    1414   !!             -   !  2005-11  (G. Madec) style & small optimisation 
    1515   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    … …  
    5151   LOGICAL , PUBLIC ::   ln_hpg_djc    = .FALSE.   !: s-coordinate (Density Jacobian with Cubic polynomial) 
    5252   LOGICAL , PUBLIC ::   ln_hpg_rot    = .FALSE.   !: s-coordinate (ROTated axes scheme) 
    53    REAL(wp), PUBLIC ::   rn_gamma      = 0.e0      !: weighting coefficient 
     53   REAL(wp), PUBLIC ::   rn_gamma      = 0._wp     !: weighting coefficient 
    5454   LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag 
    5555   INTEGER , PUBLIC ::   nn_dynhpg_rst = 0         !: add dynhpg implicit variables in restart ot not 
    … …  
    6363   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    6464   !! $Id$ 
    65    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     65   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6666   !!---------------------------------------------------------------------- 
    67  
    6867CONTAINS 
    6968 
    … …  
    8281      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D temporary workspace 
    8382      !!---------------------------------------------------------------------- 
    84     
     83      ! 
    8584      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn) 
    8685         ztrdu(:,:,:) = ua(:,:,:)   
    8786         ztrdv(:,:,:) = va(:,:,:)  
    8887      ENDIF       
    89  
     88      ! 
    9089      SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation 
    9190      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    … …  
    9796      CASE (  6 )   ;   CALL hpg_rot    ( kt )      ! s-coordinate (ROTated axes scheme) 
    9897      END SELECT 
    99  
     98      ! 
    10099      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 
    101100         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    … …  
    207206       
    208207      ! Local constant initialization  
    209       zcoef0 = - grav * 0.5 
     208      zcoef0 = - grav * 0.5_wp 
    210209 
    211210      ! Surface value 
    … …  
    270269 
    271270      ! Local constant initialization 
    272       zcoef0 = - grav * 0.5 
    273  
    274       !  Surface value 
     271      zcoef0 = - grav * 0.5_wp 
     272 
     273      !  Surface value (also valid in partial step case) 
    275274      DO jj = 2, jpjm1 
    276275         DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    305304      END DO 
    306305 
    307       ! partial steps correction at the last level  (new gradient with  intgrd.F) 
     306      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90) 
    308307# if defined key_vectopt_loop 
    309308         jj = 1 
    … …  
    313312         DO ji = 2, jpim1 
    314313# endif 
    315             iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1 
    316             ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1 
     314            iku = mbku(ji,jj) 
     315            ikv = mbkv(ji,jj) 
    317316            zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) ) 
    318317            zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) ) 
    319             ! on i-direction 
    320             IF ( iku > 2 ) THEN 
    321                ! subtract old value 
    322                ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) 
    323                ! compute the new one 
    324                zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1)   & 
    325                   + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 
    326                ! add the new one to the general momentum trend 
    327                ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) 
     318            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more) 
     319               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value 
     320               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one 
     321                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 
     322               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
    328323            ENDIF 
    329             ! on j-direction 
    330             IF ( ikv > 2 ) THEN 
    331                ! subtract old value 
    332                va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) 
    333                ! compute the new one 
    334                zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1)   & 
    335                   + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 
    336                ! add the new one to the general momentum trend 
    337                va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) 
     324            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more) 
     325               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value 
     326               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one 
     327                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 
     328               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
    338329            ENDIF 
    339330# if ! defined key_vectopt_loop 
    … …  
    379370 
    380371      ! Local constant initialization 
    381       zcoef0 = - grav * 0.5 
     372      zcoef0 = - grav * 0.5_wp 
    382373      ! To use density and not density anomaly 
    383       IF ( lk_vvl ) THEN   ;     znad = 1.            ! Variable volume 
    384       ELSE                 ;     znad = 0.e0          ! Fixed volume 
     374      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
     375      ELSE                 ;     znad = 0._wp         ! Fixed volume 
    385376      ENDIF 
    386377 
    … …  
    394385               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
    395386            ! s-coordinate pressure gradient correction 
    396             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2*znad )   & 
     387            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    397388               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    398             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2*znad )   & 
     389            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    399390               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    400391            ! add to the general momentum trend 
    … …  
    416407                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    417408               ! s-coordinate pressure gradient correction 
    418                zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2*znad )   & 
     409               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    419410                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
    420                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2*znad )   & 
     411               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    421412                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
    422413               ! add to the general momentum trend 
    … …  
    465456 
    466457      ! Local constant initialization 
    467       zcoef0 = - grav * 0.5 
     458      zcoef0 = - grav * 0.5_wp 
    468459  
    469460      ! Surface value 
    … …  
    497488                  &                                     -fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1) ) 
    498489               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 
    499                   &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk)   & 
     490                  &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk)     & 
    500491                  &                                     -fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk)   ) & 
    501                   &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) & 
     492                  &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   & 
    502493                  &                                     -fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1) ) 
    503494               ! s-coordinate pressure gradient correction 
    … …  
    543534 
    544535      ! Local constant initialization 
    545       zcoef0 = - grav * 0.5 
    546       zalph  = 0.5 - rn_gamma    ! weighting coefficients (alpha=0.5-rn_gamma 
    547       zbeta  = 0.5 + rn_gamma    !                        (beta =1-alpha=0.5+rn_gamma 
     536      zcoef0 = - grav * 0.5_wp 
     537      zalph  = 0.5_wp - rn_gamma    ! weighting coefficients (alpha=0.5-rn_gamma 
     538      zbeta  = 0.5_wp + rn_gamma    !                        (beta =1-alpha=0.5+rn_gamma 
    548539 
    549540      ! Surface value (no ponderation) 
    … …  
    628619 
    629620      ! Local constant initialization 
    630       zcoef0 = - grav * 0.5 
    631       z1_10  = 1.0 / 10.0 
    632       z1_12  = 1.0 / 12.0 
     621      zcoef0 = - grav * 0.5_wp 
     622      z1_10  = 1._wp / 10._wp 
     623      z1_12  = 1._wp / 12._wp 
    633624 
    634625      !---------------------------------------------------------------------------------------- 
    … …  
    662653         DO jj = 2, jpjm1 
    663654            DO ji = fs_2, fs_jpim1   ! vector opt. 
    664                cffw = 2.0 * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1) 
    665  
    666                cffu = 2.0 * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  ) 
    667                cffx = 2.0 * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  ) 
     655               cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1) 
     656 
     657               cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  ) 
     658               cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  ) 
    668659   
    669                cffv = 2.0 * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  ) 
    670                cffy = 2.0 * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  ) 
     660               cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  ) 
     661               cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  ) 
    671662 
    672663               IF( cffw > zep) THEN 
    673                   drhow(ji,jj,jk) = 2.0 *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   & 
    674                      &                  / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 
     664                  drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   & 
     665                     &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 
    675666               ELSE 
    676                   drhow(ji,jj,jk) = 0.e0 
     667                  drhow(ji,jj,jk) = 0._wp 
    677668               ENDIF 
    678669 
    679                dzw(ji,jj,jk) = 2.0 *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   & 
    680                   &                / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 
     670               dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   & 
     671                  &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 
    681672 
    682673               IF( cffu > zep ) THEN 
    683                   drhou(ji,jj,jk) = 2.0 *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   & 
    684                      &                  / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 
     674                  drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   & 
     675                     &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 
    685676               ELSE 
    686                   drhou(ji,jj,jk ) = 0.e0 
     677                  drhou(ji,jj,jk ) = 0._wp 
    687678               ENDIF 
    688679 
    689680               IF( cffx > zep ) THEN 
    690                   dzu(ji,jj,jk) = 2.0*dzx(ji+1,jj,jk)*dzx(ji,jj,jk)   & 
    691                      &            /(dzx(ji+1,jj,jk)+dzx(ji,jj,jk)) 
     681                  dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   & 
     682                     &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 
    692683               ELSE 
    693                   dzu(ji,jj,jk) = 0.e0 
     684                  dzu(ji,jj,jk) = 0._wp 
    694685               ENDIF 
    695686 
    696687               IF( cffv > zep ) THEN 
    697                   drhov(ji,jj,jk) = 2.0 *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   & 
    698                      &                  / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 
     688                  drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   & 
     689                     &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 
    699690               ELSE 
    700                   drhov(ji,jj,jk) = 0.e0 
     691                  drhov(ji,jj,jk) = 0._wp 
    701692               ENDIF 
    702693 
    703694               IF( cffy > zep ) THEN 
    704                   dzv(ji,jj,jk) = 2.0 *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   & 
    705                      &                / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 
     695                  dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   & 
     696                     &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 
    706697               ELSE 
    707                   dzv(ji,jj,jk) = 0.e0 
     698                  dzv(ji,jj,jk) = 0._wp 
    708699               ENDIF 
    709700 
    … …  
    715706      ! apply boundary conditions at top and bottom using 5.36-5.37 
    716707      !---------------------------------------------------------------------------------- 
    717       drhow(:,:, 1 ) = 1.5 * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5 * drhow(:,:,  2  ) 
    718       drhou(:,:, 1 ) = 1.5 * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5 * drhou(:,:,  2  ) 
    719       drhov(:,:, 1 ) = 1.5 * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5 * drhov(:,:,  2  ) 
    720  
    721       drhow(:,:,jpk) = 1.5 * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5 * drhow(:,:,jpkm1) 
    722       drhou(:,:,jpk) = 1.5 * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5 * drhou(:,:,jpkm1) 
    723       drhov(:,:,jpk) = 1.5 * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5 * drhov(:,:,jpkm1) 
     708      drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  ) 
     709      drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  ) 
     710      drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  ) 
     711 
     712      drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 
     713      drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 
     714      drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 
    724715 
    725716 
    … …  
    733724      DO jj = 2, jpjm1 
    734725         DO ji = fs_2, fs_jpim1   ! vector opt. 
    735             rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )            & 
    736                &                   * (  rhd(ji,jj,1)                                 & 
    737                &                     + 0.5 * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         & 
    738                &                           * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   & 
    739                &                           / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  )  
     726            rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )               & 
     727               &                   * (  rhd(ji,jj,1)                                    & 
     728               &                     + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         & 
     729               &                              * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   & 
     730               &                              / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  )  
    740731         END DO 
    741732      END DO 
    … …  
    849840      !  Local constant initialization 
    850841      ! ------------------------------- 
    851       zcoef0 = - grav * 0.5 
    852       zforg  = 0.95e0 
    853       zfrot  = 1.e0 - zforg 
     842      zcoef0 = - grav * 0.5_wp 
     843      zforg  = 0.95_wp 
     844      zfrot  = 1._wp - zforg 
    854845 
    855846      ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r2287 r2450  
    6969      REAL(wp) ::   zzwi, zzws, zrhs       ! temporary scalars 
    7070      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! 3D workspace 
    71 #if defined key_zdfgls 
    72       INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 
    73       REAL(wp) :: zcbcu, zcbcv 
    74 #endif 
    7571      !!---------------------------------------------------------------------- 
    7672 
    … …  
    168164      DO jj = 2, jpjm1 
    169165         DO ji = fs_2, fs_jpim1   ! vector opt. 
    170             ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    171             ikbum1 = MAX( ikbu-1, 1 ) 
    172             wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,ikbum1) * umask(ji,jj,ikbum1) 
     166            wbotu(ji,jj) = bfrua(ji,jj) * ua(ji,jj,mbku(ji,jj)) * umask(ji,jj,mbku(ji,jj)) 
    173167         END DO 
    174168      END DO 
    … …  
    269263      DO jj = 2, jpjm1 
    270264         DO ji = fs_2, fs_jpim1   ! vector opt. 
    271             ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    272             ikbvm1 = MAX( ikbv-1, 1 ) 
    273             wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1) 
     265            wbotv(ji,jj) = bfrva(ji,jj) * va(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,mbku(ji,jj)) 
    274266         END DO 
    275267      END DO 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r2371 r2450  
    44   !! Ocean physics:  variable eddy induced velocity coefficients 
    55   !!====================================================================== 
     6   !! History :  OPA  ! 1999-03  (G. Madec, A. Jouzeau)  Original code 
     7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  Free form, F90 
     8   !!---------------------------------------------------------------------- 
    69#if   defined key_traldf_eiv   &&   defined key_traldf_c2d 
    710   !!---------------------------------------------------------------------- 
    … …  
    1114   !!   ldf_eiv      : compute the eddy induced velocity coefficients 
    1215   !!---------------------------------------------------------------------- 
    13    !! * Modules used 
    1416   USE oce             ! ocean dynamics and tracers 
    1517   USE dom_oce         ! ocean space and time domain 
    … …  
    2224   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2325   USE prtctl          ! Print control 
    24    USE iom 
     26   USE iom             ! I/O library 
    2527 
    2628   IMPLICIT NONE 
    2729   PRIVATE 
    2830    
    29    !! * Routine accessibility 
    30    PUBLIC ldf_eiv               ! routine called by step.F90 
    31    !!---------------------------------------------------------------------- 
    32    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    33    !! $Id$ 
    34    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    35    !!---------------------------------------------------------------------- 
     31   PUBLIC   ldf_eiv    ! routine called by step.F90 
     32    
    3633   !! * Substitutions 
    3734#  include "domzgr_substitute.h90" 
    3835#  include "vectopt_loop_substitute.h90" 
    3936   !!---------------------------------------------------------------------- 
    40  
     37   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     38   !! $Id$ 
     39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     40   !!---------------------------------------------------------------------- 
    4141CONTAINS 
    4242 
    … …  
    4646      !! 
    4747      !! ** Purpose :   Compute the eddy induced velocity coefficient from the 
    48       !!      growth rate of baroclinic instability. 
     48      !!              growth rate of baroclinic instability. 
    4949      !! 
    5050      !! ** Method : 
    5151      !! 
    52       !! ** Action : - uslp(),  : i- and j-slopes of neutral surfaces 
    53       !!             - vslp()      at u- and v-points, resp. 
    54       !!             - wslpi(),  : i- and j-slopes of neutral surfaces 
    55       !!             - wslpj()     at w-points.  
    56       !! 
    57       !! History : 
    58       !!   8.1  !  99-03  (G. Madec, A. Jouzeau)  Original code 
    59       !!   8.5  !  02-06  (G. Madec)  Free form, F90 
     52      !! ** Action : - uslp , vslp  : i- and j-slopes of neutral surfaces at u- & v-points 
     53      !!             - wslpi, wslpj : i- and j-slopes of neutral surfaces at w-points.  
    6054      !!---------------------------------------------------------------------- 
    61       !! * Arguments 
    62       INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx 
    63        
    64       !! * Local declarations 
    65       INTEGER ::   ji, jj, jk           ! dummy loop indices 
    66       REAL(wp) ::   & 
    67          zfw, ze3w, zn2, zf20,       &  ! temporary scalars 
    68          zaht, zaht_min 
    69       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    70          zn, zah, zhw, zross            ! workspace 
     55      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx 
     56      !! 
     57      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     58      REAL(wp) ::   zfw, ze3w, zn2, zf20, zaht, zaht_min      ! temporary scalars 
     59      REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zross   ! 2D workspace 
    7160      !!---------------------------------------------------------------------- 
    7261       
    … …  
    7968      ! 0. Local initialization 
    8069      ! ----------------------- 
    81       zn   (:,:) = 0.e0 
    82       zhw  (:,:) = 5.e0 
    83       zah  (:,:) = 0.e0 
    84       zross(:,:) = 0.e0 
     70      zn   (:,:) = 0._wp 
     71      zhw  (:,:) = 5._wp 
     72      zah  (:,:) = 0._wp 
     73      zross(:,:) = 0._wp 
    8574 
    8675 
    8776      ! 1. Compute lateral diffusive coefficient  
    8877      ! ---------------------------------------- 
    89       IF (ln_traldf_grif) THEN 
     78      IF( ln_traldf_grif ) THEN 
    9079         DO jk = 1, jpk 
    9180#  if defined key_vectopt_loop   
    92             !CDIR NOVERRCHK  
     81!CDIR NOVERRCHK  
    9382            DO ji = 1, jpij   ! vector opt. 
    9483               ! Take the max of N^2 and zero then take the vertical sum 
    9584               ! of the square root of the resulting N^2 ( required to compute 
    9685               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
    97                zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 
     86               zn2 = MAX( rn2b(ji,1,jk), 0._wp ) 
    9887               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 
    9988               ! Compute elements required for the inverse time scale of baroclinic 
    … …  
    10695#  else 
    10796            DO jj = 2, jpjm1 
    108                !CDIR NOVERRCHK  
     97!CDIR NOVERRCHK  
    10998               DO ji = 2, jpim1 
    11099                  ! Take the max of N^2 and zero then take the vertical sum  
    111100                  ! of the square root of the resulting N^2 ( required to compute  
    112101                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
    113                   zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 
     102                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 
    114103                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 
    115104                  ! Compute elements required for the inverse time scale of baroclinic 
    … …  
    126115         DO jk = 1, jpk 
    127116#  if defined key_vectopt_loop   
    128             !CDIR NOVERRCHK  
     117!CDIR NOVERRCHK  
    129118            DO ji = 1, jpij   ! vector opt. 
    130119               ! Take the max of N^2 and zero then take the vertical sum 
    131120               ! of the square root of the resulting N^2 ( required to compute 
    132121               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
    133                zn2 = MAX( rn2b(ji,1,jk), 0.e0 ) 
     122               zn2 = MAX( rn2b(ji,1,jk), 0._wp ) 
    134123               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk) 
    135124               ! Compute elements required for the inverse time scale of baroclinic 
    … …  
    137126               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    138127               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk) 
    139                zah(ji,1) = zah(ji,1) + zn2   & 
    140                     * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    & 
    141                     + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   & 
    142                     * ze3w 
     128               zah(ji,1) = zah(ji,1) + zn2 * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)   & 
     129                  &                          + wslpj(ji,1,jk) * wslpj(ji,1,jk) ) * ze3w 
    143130               zhw(ji,1) = zhw(ji,1) + ze3w 
    144131            END DO 
    145132#  else 
    146133            DO jj = 2, jpjm1 
    147                !CDIR NOVERRCHK  
     134!CDIR NOVERRCHK  
    148135               DO ji = 2, jpim1 
    149136                  ! Take the max of N^2 and zero then take the vertical sum  
    150137                  ! of the square root of the resulting N^2 ( required to compute  
    151138                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f  
    152                   zn2 = MAX( rn2b(ji,jj,jk), 0.e0 ) 
     139                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 
    153140                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk) 
    154141                  ! Compute elements required for the inverse time scale of baroclinic 
    … …  
    156143                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    157144                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
    158                   zah(ji,jj) = zah(ji,jj) + zn2   & 
    159                        * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    & 
    160                        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  & 
    161                        * ze3w 
     145                  zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     146                     &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 
    162147                  zhw(ji,jj) = zhw(ji,jj) + ze3w 
    163148               END DO 
    … …  
    178163      END DO 
    179164 
    180       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R02 
     165      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2 
    181166         DO jj = 2, jpjm1 
    182167!CDIR NOVERRCHK  
    183168            DO ji = fs_2, fs_jpim1   ! vector opt. 
    184                ! Take the minimum between aeiw and aeiv0 for depth levels 
    185                ! lower than 20 (21 in w- point) 
    186                IF( mbathy(ji,jj) <= 21. ) aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. ) 
     169               ! Take the minimum between aeiw and 1000 m2/s over shelves (depth shallower than 650 m) 
     170               IF( mbkt(ji,jj) <= 20 )   aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. ) 
    187171            END DO 
    188172         END DO 
    … …  
    190174 
    191175      ! Decrease the coefficient in the tropics (20N-20S)  
    192       zf20 = 2. * omega * sin( rad * 20. ) 
     176      zf20 = 2._wp * omega * sin( rad * 20._wp ) 
    193177      DO jj = 2, jpjm1 
    194178         DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    205189         END DO 
    206190      ENDIF 
    207  
    208       ! lateral boundary condition on aeiw  
    209       CALL lbc_lnk( aeiw, 'W', 1. ) 
     191      CALL lbc_lnk( aeiw, 'W', 1. )      ! lateral boundary condition on aeiw  
     192 
    210193 
    211194      ! Average the diffusive coefficient at u- v- points  
    212195      DO jj = 2, jpjm1 
    213196         DO ji = fs_2, fs_jpim1   ! vector opt. 
    214             aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) ) 
    215             aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) ) 
     197            aeiu(ji,jj) = 0.5_wp * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) ) 
     198            aeiv(ji,jj) = 0.5_wp * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) ) 
    216199         END DO  
    217200      END DO  
    218  
    219       ! lateral boundary condition on aeiu, aeiv 
    220       CALL lbc_lnk( aeiu, 'U', 1. ) 
    221       CALL lbc_lnk( aeiv, 'V', 1. ) 
    222  
    223       IF(ln_ctl)   THEN 
     201      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )      ! lateral boundary condition 
     202 
     203 
     204      IF(ln_ctl) THEN 
    224205         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1) 
    225206         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1) 
    … …  
    228209      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth) 
    229210      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN 
    230          zf20     = 2. * omega * SIN( rad * 20. ) 
    231          zaht_min = 100.                              ! minimum value for aht 
     211         zf20     = 2._wp * omega * SIN( rad * 20._wp ) 
     212         zaht_min = 100._wp                           ! minimum value for aht 
    232213         DO jj = 1, jpj 
    233214            DO ji = 1, jpi 
    234                zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  & 
     215               zaht      = ( 1._wp -  MIN( 1._wp , ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  & 
    235216                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths 
    236217               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 ) 
    … …  
    246227      ENDIF 
    247228 
    248       IF( aeiv0 == 0.e0 ) THEN 
    249          aeiu(:,:) = 0.e0 
    250          aeiv(:,:) = 0.e0 
    251          aeiw(:,:) = 0.e0 
     229      IF( aeiv0 == 0._wp ) THEN 
     230         aeiu(:,:) = 0._wp 
     231         aeiv(:,:) = 0._wp 
     232         aeiw(:,:) = 0._wp 
    252233      ENDIF 
    253234 
    254235      CALL iom_put( "aht2d"    , ahtw )   ! lateral eddy diffusivity 
    255236      CALL iom_put( "aht2d_eiv", aeiw )   ! EIV lateral eddy diffusivity 
    256  
     237      ! 
    257238   END SUBROUTINE ldf_eiv 
    258239 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r2435 r2450  
    77   !!            8.0  ! 1997-06  (G. Madec)  optimization, lbc 
    88   !!            8.1  ! 1999-10  (A. Jouzeau)  NEW profile in the mixed layer 
    9    !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90 
    10    !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates 
     9   !!   NEMO     1.0  ! 2002-10  (G. Madec)  Free form, F90 
     10   !!             -   ! 2005-10  (A. Beckmann)  correction for s-coordinates 
    1111   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator 
    1212   !!             -   ! 2010-11  (F. Dupond, G. Madec)  bug correction in slopes just below the ML 
    … …  
    116116      zm1_2g = -0.5_wp / grav 
    117117      ! 
    118       zww(:,:,:) = 0.e0 
    119       zwz(:,:,:) = 0.e0 
     118      zww(:,:,:) = 0._wp 
     119      zwz(:,:,:) = 0._wp 
    120120      ! 
    121121      DO jk = 1, jpk             !==   i- & j-gradient of density   ==! 
    … …  
    135135            DO ji = 1, jpim1 
    136136# endif 
    137                iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1               ! last ocean level 
    138                ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1 
    139                zgru(ji,jj,iku) = gru(ji,jj)  
    140                zgrv(ji,jj,ikv) = grv(ji,jj)                
     137               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)  
     138               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)                
    141139            END DO 
    142140         END DO 
    … …  
    329327         !                                                    ! Gibraltar Strait 
    330328         ij0 =  50   ;   ij1 =  53 
    331          ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
     329         ii0 =  69   ;   ii1 =  71   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    332330         ij0 =  51   ;   ij1 =  53 
    333          ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
    334          ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
    335          ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
     331         ii0 =  68   ;   ii1 =  71   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
     332         ii0 =  69   ;   ii1 =  71   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
     333         ii0 =  69   ;   ii1 =  71   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    336334         ! 
    337335         !                                                    ! Mediterrannean Sea 
    338336         ij0 =  49   ;   ij1 =  56 
    339          ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
     337         ii0 =  71   ;   ii1 =  90   ;   uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    340338         ij0 =  50   ;   ij1 =  56 
    341          ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
    342          ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
    343          ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0.e0 
     339         ii0 =  70   ;   ii1 =  90   ;   vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
     340         ii0 =  71   ;   ii1 =  90   ;   wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
     341         ii0 =  71   ;   ii1 =  90   ;   wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 
    344342      ENDIF 
    345343 
    … …  
    403401         DO jj = 1, jpjm1 
    404402            DO ji = 1, fs_jpim1   ! vector opt. 
    405                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 
     403               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 
    406404               zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 
    407                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 
     405               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 
    408406               zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 
    409407            END DO 
    … …  
    418416            DO ji = 1, jpim1 
    419417# endif 
    420                iku = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1, 2 )    ! last ocean level 
    421                ikv = MAX( MIN( mbathy(ji,jj), mbathy(ji,jj+1  ) ) - 1, 2 ) 
    422                zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem)                          ! i-gradient of T and S 
    423                zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal) 
    424                zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem)                          ! j-gradient of T and S 
    425                zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal) 
     418               zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem)                           ! i-gradient of T and S 
     419               zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal) 
     420               zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem)                           ! j-gradient of T and S 
     421               zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal) 
    426422            END DO 
    427423         END DO 
    … …  
    453449            DO jj = 1, jpj                       ! NB: not masked due to the minimum value set 
    454450               DO ji = 1, jpi   ! vector opt.  
    455                    zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) )   & 
     451                  zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) )   & 
    456452                     &       / fse3w(ji,jj,jk+kp) 
    457453                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )                    ! force zdzrho >= repsln 
    … …  
    462458      ! 
    463459      DO jj = 1, jpj                         !==  Reciprocal depth of the w-point below ML base  ==! 
    464          DO ji = 1, jpi                            ! i.e. 1 / (hmld+e3t(nmln))  where hmld=depw(nmln) 
    465             jk = MIN( nmln(ji,jj), mbathy(ji,jj) - 1 ) + 1     ! MIN in case ML depth is the ocean depth 
     460         DO ji = 1, jpi 
     461            jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth 
    466462            z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk) 
    467463         END DO 
    … …  
    471467      ! 
    472468      wslp2  (:,:,:)     = 0._wp                                           ! wslp2 will be cumulated 3D field set to zero 
    473       triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero 
     469      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp      ! set surface and bottom slope to zero 
    474470      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp 
    475471!!gm _iso set to zero missing 
    … …  
    486482               DO ji = 1, fs_jpim1 
    487483                  ip = jl   ;   jp = jl 
    488                   jk = MIN( nmln(ji+ip,jj), mbathy(ji+ip,jj) - 1 ) + 1     ! ML level+1 (MIN in case ML depth is the ocean depth) 
     484                  jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1         ! ML level+1 (MIN in case ML depth is the ocean depth) 
    489485                  zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      & 
    490486                     &      + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk) 
    491                   jk = MIN( nmln(ji,jj+jp), mbathy(ji,jj+jp) - 1 ) + 1 
     487                  jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 
    492488                  ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      & 
    493489                     &      + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk) 
    … …  
    612608      zm1_2g = -0.5_wp / grav 
    613609      ! 
    614       uslpml (1,:) = 0.e0      ;      uslpml (jpi,:) = 0.e0 
    615       vslpml (1,:) = 0.e0      ;      vslpml (jpi,:) = 0.e0 
    616       wslpiml(1,:) = 0.e0      ;      wslpiml(jpi,:) = 0.e0 
    617       wslpjml(1,:) = 0.e0      ;      wslpjml(jpi,:) = 0.e0 
     610      uslpml (1,:) = 0._wp      ;      uslpml (jpi,:) = 0._wp 
     611      vslpml (1,:) = 0._wp      ;      vslpml (jpi,:) = 0._wp 
     612      wslpiml(1,:) = 0._wp      ;      wslpiml(jpi,:) = 0._wp 
     613      wslpjml(1,:) = 0._wp      ;      wslpjml(jpi,:) = 0._wp 
    618614      ! 
    619615      !                          !==   surface mixed layer mask   ! 
    … …  
    627623# endif 
    628624               ik = nmln(ji,jj) - 1 
    629                IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1.e0 
    630                ELSE                  ;   omlmask(ji,jj,jk) = 0.e0 
     625               IF( jk <= ik ) THEN   ;   omlmask(ji,jj,jk) = 1._wp 
     626               ELSE                  ;   omlmask(ji,jj,jk) = 0._wp 
    631627               ENDIF 
    632628            END DO 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r2384 r2450  
    1919   USE phycst          ! physical constants 
    2020   USE sbc_oce         ! surface boundary condition variables 
    21    USE fldread         ! ??? 
     21   USE fldread         ! read input field at current time step 
    2222   USE in_out_manager  ! I/O manager 
    2323   USE iom             ! I/O module 
    2424   USE restart         ! restart 
    25    USE closea 
     25   USE closea          ! closed seas 
    2626 
    2727   IMPLICIT NONE 
    2828   PRIVATE 
    2929 
    30    PUBLIC sbc_rnf          ! routine call in sbcmod module 
    31    PUBLIC sbc_rnf_div      ! routine called in sshwzv module 
     30   PUBLIC   sbc_rnf       ! routine call in sbcmod module 
     31   PUBLIC   sbc_rnf_div   ! routine called in sshwzv module 
    3232 
    3333   !                                                      !!* namsbc_rnf namelist * 
    … …  
    4343   TYPE(FLD_N)                ::   sn_dep_rnf             !: information about the depth which river inflow affects 
    4444   LOGICAL           , PUBLIC ::   ln_rnf_mouth = .false. !: specific treatment in mouths vicinity 
    45    REAL(wp)          , PUBLIC ::   rn_hrnf      = 0.e0    !: runoffs, depth over which enhanced vertical mixing is used 
    46    REAL(wp)          , PUBLIC ::   rn_avt_rnf   = 0.e0    !: runoffs, value of the additional vertical mixing coef. [m2/s] 
    47    REAL(wp)          , PUBLIC ::   rn_rfact     = 1.e0    !: multiplicative factor for runoff 
    48  
    49    INTEGER , PUBLIC                     ::   nkrnf = 0    !: number of levels over which Kz is increased at river mouths 
    50    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rnfmsk       !: river mouth mask (hori.) 
    51    REAL(wp), PUBLIC, DIMENSION(jpk)     ::   rnfmsk_z     !: river mouth mask (vert.) 
    52  
    53    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf       !: structure: river runoff (file information, fields read) 
    54  
    55    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     !: structure: river runoff salinity (file information, fields read)   
    56    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     !: structure: river runoff temperature (file information, fields read)   
     45   REAL(wp)          , PUBLIC ::   rn_hrnf      = 0._wp   !: runoffs, depth over which enhanced vertical mixing is used 
     46   REAL(wp)          , PUBLIC ::   rn_avt_rnf   = 0._wp   !: runoffs, value of the additional vertical mixing coef. [m2/s] 
     47   REAL(wp)          , PUBLIC ::   rn_rfact     = 1._wp   !: multiplicative factor for runoff 
     48 
     49   INTEGER , PUBLIC                          ::   nkrnf = 0         !: nb of levels over which Kz is increased at river mouths 
     50   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)      ::   rnfmsk            !: river mouth mask (hori.) 
     51   REAL(wp), PUBLIC, DIMENSION(jpk)          ::   rnfmsk_z          !: river mouth mask (vert.) 
     52   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)      ::   h_rnf             !: depth of runoff in m 
     53   INTEGER,  PUBLIC, DIMENSION(jpi,jpj)      ::   nk_rnf            !: depth of runoff in model levels 
     54   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc  !: before and now T & S contents of runoffs  [K.m/s & PSU.m/s] 
     55    
     56   REAL(wp) ::   r1_rau0   ! = 1 / rau0  
     57 
     58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf       ! structure: river runoff (file information, fields read) 
     59   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     ! structure: river runoff salinity (file information, fields read)   
     60   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     ! structure: river runoff temperature (file information, fields read)   
    5761  
    58    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   h_rnf        !: depth of runoff in m 
    59    INTEGER,  PUBLIC, DIMENSION(jpi,jpj) ::   nk_rnf       !: depth of runoff in model levels 
    60  
    61    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc  !: before and now T & S contents of runoffs  [K.m/s & PSU.m/s] 
    62  
    6362   !! * Substitutions   
    6463#  include "domzgr_substitute.h90"   
    … …  
    8584      !! 
    8685      INTEGER  ::   ji, jj   ! dummy loop indices 
    87       REAL(wp) ::   z1_rau0  ! local scalar 
    8886      !!---------------------------------------------------------------------- 
    8987      !                                    
    … …  
    105103         IF( ln_rnf_tem  )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required 
    106104         IF( ln_rnf_sal  )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required 
    107  
     105         ! 
    108106         ! Runoff reduction only associated to the ORCA2_LIM configuration 
    109107         ! when reading the NetCDF file runoff_1m_nomask.nc 
    110108         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 )   THEN 
    111             DO jj = 1, jpj 
    112                DO ji = 1, jpi 
    113                   IF( gphit(ji,jj) > 40 .AND. gphit(ji,jj) < 65 )   sf_rnf(1)%fnow(ji,jj,1) = 0.85 * sf_rnf(1)%fnow(ji,jj,1) 
    114                END DO 
    115             END DO 
    116          ENDIF 
    117  
    118          ! C a u t i o n : runoff is negative and in kg/m2/s  
     109            WHERE( 40._wp < gphit(:,:) .AND. gphit(:,:) < 65._wp ) 
     110               sf_rnf(1)%fnow(:,:,1) = 0.85 * sf_rnf(1)%fnow(:,:,1) 
     111            END WHERE 
     112         ENDIF 
     113         ! 
    119114         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    120             rnf(:,:)  = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )   
     115            rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )   
    121116            ! 
    122             z1_rau0 = 1.e0 / rau0 
    123             !                                                              ! set temperature & salinity content of runoffs 
    124             IF( ln_rnf_tem )   THEN                                        ! use runoffs temperature data 
    125                rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0 
    126                WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 )                      ! if missing data value use SST as runoffs temperature   
    127                    rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
    128                ENDWHERE 
    129             ELSE                                                           ! use SST as runoffs temperature 
    130                rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
     117            r1_rau0 = 1._wp / rau0 
     118            !                                                     ! set temperature & salinity content of runoffs 
     119            IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data 
     120               rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 
     121               WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999 )                 ! if missing data value use SST as runoffs temperature   
     122                   rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
     123               END WHERE 
     124            ELSE                                                        ! use SST as runoffs temperature 
     125               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
    131126            ENDIF   
    132             !                                                              ! use runoffs salinity data  
    133             IF( ln_rnf_sal ) rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * z1_rau0 
    134             !                                                              ! else use S=0 for runoffs (done one for all in the init) 
     127            !                                                           ! use runoffs salinity data  
     128            IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 
     129            !                                                           ! else use S=0 for runoffs (done one for all in the init) 
    135130            ! 
    136             IF( ln_rnf_tem .OR. ln_rnf_sal ) THEN                         ! runoffs as outflow: use ocean SST and SSS 
    137                WHERE( rnf(:,:) < 0.e0 )                                    ! example baltic model when flow is out of domain  
    138                   rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * z1_rau0 
    139                   rnf_tsc(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * z1_rau0 
    140                ENDWHERE 
     131            IF( ln_rnf_tem .OR. ln_rnf_sal ) THEN                 ! runoffs as outflow: use ocean SST and SSS 
     132               WHERE( rnf(:,:) < 0._wp )                                 ! example baltic model when flow is out of domain  
     133                  rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
     134                  rnf_tsc(:,:,jp_sal) = sss_m(:,:) * rnf(:,:) * r1_rau0 
     135               END WHERE 
    141136            ENDIF 
    142137            ! 
    … …  
    190185      !! 
    191186      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    192       REAL(wp) ::   z1_rau0   ! local scalar 
     187      REAL(wp) ::   r1_rau0   ! local scalar 
    193188      REAL(wp) ::   zfact     ! local scalar 
    194189      !!---------------------------------------------------------------------- 
    195190      ! 
    196       zfact = 0.5e0 
    197       ! 
    198       z1_rau0 = 1.e0 / rau0 
     191      zfact = 0.5_wp 
     192      ! 
     193      r1_rau0 = 1._wp / rau0 
    199194      IF( ln_rnf_depth ) THEN      !==   runoff distributed over several levels   ==! 
    200195         IF( lk_vvl ) THEN             ! variable volume case  
    201196            DO jj = 1, jpj                   ! update the depth over which runoffs are distributed 
    202197               DO ji = 1, jpi 
    203                   h_rnf(ji,jj) = 0.e0  
     198                  h_rnf(ji,jj) = 0._wp  
    204199                  DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
    205200                     h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box  
    … …  
    207202                  !                          ! apply the runoff input flow 
    208203                  DO jk = 1, nk_rnf(ji,jj) 
    209                      phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * z1_rau0 / h_rnf(ji,jj) 
     204                     phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj) 
    210205                  END DO 
    211206               END DO 
    … …  
    215210               DO ji = 1, jpi 
    216211                  DO jk = 1, nk_rnf(ji,jj) 
    217                      phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * z1_rau0 / h_rnf(ji,jj) 
     212                     phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rau0 / h_rnf(ji,jj) 
    218213                  END DO 
    219214               END DO 
    … …  
    224219            h_rnf(:,:) = fse3t(:,:,1)   ! recalculate h_rnf to be depth of top box 
    225220         ENDIF 
    226          phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * z1_rau0 / fse3t(:,:,1) 
     221         phdivn(:,:,1) = phdivn(:,:,1) - ( rnf(:,:) + rnf_b(:,:) ) * zfact * r1_rau0 / fse3t(:,:,1) 
    227222      ENDIF 
    228223      ! 
    … …  
    303298         CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' ) 
    304299         ! 
    305          IF( ln_rnf_tem ) THEN                     ! Create (if required) sf_t_rnf structure 
     300         IF( ln_rnf_tem ) THEN                      ! Create (if required) sf_t_rnf structure 
    306301            IF(lwp) WRITE(numout,*) 
    307302            IF(lwp) WRITE(numout,*) '          runoffs temperatures read in a file' 
    … …  
    326321            CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )   
    327322         ENDIF 
    328  
    329   
    330          IF ( ln_rnf_depth ) THEN                     ! depth of runoffs set from a file  
     323         ! 
     324         IF( ln_rnf_depth ) THEN                    ! depth of runoffs set from a file  
    331325            IF(lwp) WRITE(numout,*) 
    332326            IF(lwp) WRITE(numout,*) '          runoffs depth read in a file' 
    333327            rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname )   
    334328            CALL iom_open ( rn_dep_file, inum )                           ! open file   
    335             CALL iom_get  ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf )    ! read the river mouth array   
    336             CALL iom_close( inum )                                      ! close file   
    337    
    338             nk_rnf(:,:) = 0                              ! set the number of level over which river runoffs are applied 
     329            CALL iom_get  ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array   
     330            CALL iom_close( inum )                                        ! close file   
     331            ! 
     332            nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied 
    339333            DO jj = 1, jpj   
    340               DO ji = 1, jpi   
    341                 IF ( h_rnf(ji,jj) > 0.e0 ) THEN   
    342                   jk = 2   
    343                   DO WHILE ( jk /= ( mbathy(ji,jj) - 1 ) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 ;  ENDDO   
    344                   nk_rnf(ji,jj) = jk   
    345                 ELSE IF ( h_rnf(ji,jj) == -1   ) THEN   ;  nk_rnf(ji,jj) = 1   
    346                 ELSE IF ( h_rnf(ji,jj) == -999 ) THEN   ;  nk_rnf(ji,jj) = mbathy(ji,jj) - 1 
    347                 ELSE IF ( h_rnf(ji,jj) /= 0 ) THEN   
    348                   CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )   
    349                   WRITE(999,*) 'ji, jj, rnf(ji,jj) :', ji, jj, rnf(ji,jj)   
    350                 ENDIF   
    351               ENDDO   
    352             ENDDO   
    353             DO jj = 1, jpj                               ! set the associated depth  
    354               DO ji = 1, jpi  
    355                 h_rnf(ji,jj) = 0.e0 
    356                 DO jk = 1, nk_rnf(ji,jj)                         
    357                    h_rnf(ji,jj) = h_rnf(ji,jj)+fse3t(ji,jj,jk)   
    358                 ENDDO 
    359               ENDDO 
    360             ENDDO 
     334               DO ji = 1, jpi   
     335                  IF( h_rnf(ji,jj) > 0._wp ) THEN   
     336                     jk = 2   
     337                     DO WHILE ( jk /= mbkt(ji,jj) .AND. fsdept(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 ;  END DO   
     338                     nk_rnf(ji,jj) = jk   
     339                  ELSEIF( h_rnf(ji,jj) == -1   ) THEN   ;  nk_rnf(ji,jj) = 1   
     340                  ELSEIF( h_rnf(ji,jj) == -999 ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj) 
     341                  ELSEIF( h_rnf(ji,jj) /=  0  ) THEN   
     342                     CALL ctl_stop( 'runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )   
     343                     WRITE(999,*) 'ji, jj, rnf(ji,jj) :', ji, jj, rnf(ji,jj)   
     344                  ENDIF   
     345               END DO   
     346            END DO   
     347            DO jj = 1, jpj                                ! set the associated depth  
     348               DO ji = 1, jpi  
     349                  h_rnf(ji,jj) = 0._wp 
     350                  DO jk = 1, nk_rnf(ji,jj)                         
     351                     h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   
     352                  END DO 
     353               END DO 
     354            END DO 
    361355         ELSE                                       ! runoffs applied at the surface  
    362356            nk_rnf(:,:) = 1   
    363             h_rnf(:,:) = fse3t(:,:,1) 
     357            h_rnf (:,:) = fse3t(:,:,1) 
    364358         ENDIF   
    365       !  
    366       ENDIF 
    367  
    368       rnf_tsc(:,:,:) = 0.e0                 ! runoffs temperature & salinty contents initilisation 
     359         !  
     360      ENDIF 
     361      ! 
     362      rnf_tsc(:,:,:) = 0._wp                    ! runoffs temperature & salinty contents initilisation 
     363      ! 
    369364      !                                   ! ======================== 
    370365      !                                   !   River mouth vicinity 
    … …  
    380375         ! 
    381376         nkrnf = 0                                  ! Number of level over which Kz increase 
    382          IF( rn_hrnf > 0.e0 ) THEN 
     377         IF( rn_hrnf > 0._wp ) THEN 
    383378            nkrnf = 2 
    384379            DO WHILE( nkrnf /= jpkm1 .AND. gdepw_0(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1   ;   END DO 
    … …  
    398393         IF(lwp) WRITE(numout,*) 
    399394         IF(lwp) WRITE(numout,*) '          No specific treatment at river mouths' 
    400          rnfmsk  (:,:) = 0.e0  
    401          rnfmsk_z(:)   = 0.e0 
     395         rnfmsk  (:,:) = 0._wp  
     396         rnfmsk_z(:)   = 0._wp 
    402397         nkrnf = 0 
    403398      ENDIF 
    … …  
    447442      IF( nclosea == 1 )    CALL clo_rnf( rnfmsk )                ! closed sea inflow set as ruver mouth 
    448443 
    449       rnfmsk_z(:)   = 0.e0                                        ! vertical structure  
     444      rnfmsk_z(:)   = 0._wp                                        ! vertical structure  
    450445      rnfmsk_z(1)   = 1.0 
    451446      rnfmsk_z(2)   = 1.0                                         ! ********** 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r2399 r2450  
    8484      !! 
    8585      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    86       INTEGER  ::  iku, ikv         ! local integers 
    8786      REAL(wp) ::  zbtr, ztra       ! local scalars 
    8887      REAL(wp), DIMENSION(jpi,jpj) ::   zeeu, zeev, zlt   ! 2D workspace 
    … …  
    116115               END DO 
    117116            END DO 
    118             IF( ln_zps ) THEN                ! set gradient at partial step level 
     117            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level) 
    119118               DO jj = 1, jpjm1 
    120119                  DO ji = 1, jpim1 
    121                      ! last level 
    122                      iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1 
    123                      ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
    124                      IF( iku == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 
    125                      IF( ikv == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 
     120                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 
     121                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 
    126122                  END DO 
    127123               END DO 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r2399 r2450  
    101101      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
    102102      !! 
    103       INTEGER  ::  ji, jj, jk,jn   ! dummy loop indices 
    104       INTEGER  ::  iku, ikv        ! temporary integer 
     103      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    105104      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    106105      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    … …  
    141140            END DO 
    142141         END DO 
    143          IF( ln_zps ) THEN      ! partial steps correction at the last level  
     142         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level  
    144143            DO jj = 1, jpjm1 
    145144               DO ji = 1, fs_jpim1   ! vector opt. 
    146                   ! last level 
    147                   iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1 
    148                   ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
    149                   zdit(ji,jj,iku) = pgu(ji,jj,jn)           
    150                   zdjt(ji,jj,ikv) = pgv(ji,jj,jn)       
     145                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     146                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)       
    151147               END DO 
    152148            END DO 
    … …  
    215211#if defined key_diaar5 
    216212         IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    217             zztmp = 0.5 * rau0 * rcp  
    218             z2d(:,:) = 0.e0  
     213            z2d(:,:) = 0._wp  
     214            zztmp = rau0 * rcp  
    219215            DO jk = 1, jpkm1 
    220216               DO jj = 2, jpjm1 
    221217                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    222                      z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk)       & 
    223                         &       * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk)  
     218                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
    224219                  END DO 
    225220               END DO 
    226221            END DO 
     222            z2d(:,:) = zztmp * z2d(:,:) 
    227223            CALL lbc_lnk( z2d, 'U', -1. ) 
    228224            CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    229             z2d(:,:) = 0.e0  
    230225            DO jk = 1, jpkm1 
    231226               DO jj = 2, jpjm1 
    232227                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    233                      z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk)       & 
    234                         &       * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk)  
     228                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
    235229                  END DO 
    236230               END DO 
    237231            END DO 
     232            z2d(:,:) = zztmp * z2d(:,:) 
    238233            CALL lbc_lnk( z2d, 'V', -1. ) 
    239234            CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90

    r2399 r2450  
    5353  SUBROUTINE tra_ldf_iso_grif( kt, cdtype, pgu, pgv,              & 
    5454       &                                   ptb, pta, kjpt, pahtb0 ) 
    55     !!---------------------------------------------------------------------- 
    56     !!                  ***  ROUTINE tra_ldf_iso_grif  *** 
    57     !! 
    58     !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
    59     !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and  
    60     !!      add it to the general trend of tracer equation. 
    61     !! 
    62     !! ** Method  :   The horizontal component of the lateral diffusive trends  
    63     !!      is provided by a 2nd order operator rotated along neural or geopo- 
    64     !!      tential surfaces to which an eddy induced advection can be added 
    65     !!      It is computed using before fields (forward in time) and isopyc- 
    66     !!      nal or geopotential slopes computed in routine ldfslp. 
    67     !! 
    68     !!      1st part :  masked horizontal derivative of T  ( di[ t ] ) 
    69     !!      ========    with partial cell update if ln_zps=T. 
    70     !! 
    71     !!      2nd part :  horizontal fluxes of the lateral mixing operator 
    72     !!      ========     
    73     !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 
    74     !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ] 
    75     !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 
    76     !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ] 
    77     !!      take the horizontal divergence of the fluxes: 
    78     !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
    79     !!      Add this trend to the general trend (ta,sa): 
    80     !!         ta = ta + difft 
    81     !! 
    82     !!      3rd part: vertical trends of the lateral mixing operator 
    83     !!      ========  (excluding the vertical flux proportional to dk[t] ) 
    84     !!      vertical fluxes associated with the rotated lateral mixing: 
    85     !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ] 
    86     !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  } 
    87     !!      take the horizontal divergence of the fluxes: 
    88     !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
    89     !!      Add this trend to the general trend (ta,sa): 
    90     !!         pta = pta + difft 
    91     !! 
    92     !! ** Action :   Update pta arrays with the before rotated diffusion 
    93     !!---------------------------------------------------------------------- 
    94     USE oce,   zftu => ua   ! use ua as workspace 
    95     USE oce,   zftv => va   ! use va as workspace 
    96     !! 
    97     INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
    98     CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    99     INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    100     REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
    101     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    102     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    103     REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
    104     !! 
    105     INTEGER  ::  ji, jj, jk,jn   ! dummy loop indices 
    106     INTEGER  ::  ip,jp,kp        ! dummy loop indices 
    107     INTEGER  ::  iku, ikv        ! temporary integer 
    108     INTEGER  ::  ierr            ! temporary integer 
    109     REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    110     REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    111     REAL(wp) ::  zcoef0, zbtr                  !   -      - 
    112     REAL(wp), DIMENSION(jpi,jpj,0:1) ::   zdkt               ! 2D+1 workspace 
    113     REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw   ! 3D workspace 
     55      !!---------------------------------------------------------------------- 
     56      !!                  ***  ROUTINE tra_ldf_iso_grif  *** 
     57      !! 
     58      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
     59      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and  
     60      !!      add it to the general trend of tracer equation. 
     61      !! 
     62      !! ** Method  :   The horizontal component of the lateral diffusive trends  
     63      !!      is provided by a 2nd order operator rotated along neural or geopo- 
     64      !!      tential surfaces to which an eddy induced advection can be added 
     65      !!      It is computed using before fields (forward in time) and isopyc- 
     66      !!      nal or geopotential slopes computed in routine ldfslp. 
     67      !! 
     68      !!      1st part :  masked horizontal derivative of T  ( di[ t ] ) 
     69      !!      ========    with partial cell update if ln_zps=T. 
     70      !! 
     71      !!      2nd part :  horizontal fluxes of the lateral mixing operator 
     72      !!      ========     
     73      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 
     74      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ] 
     75      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ] 
     76      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ] 
     77      !!      take the horizontal divergence of the fluxes: 
     78      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  } 
     79      !!      Add this trend to the general trend (ta,sa): 
     80      !!         ta = ta + difft 
     81      !! 
     82      !!      3rd part: vertical trends of the lateral mixing operator 
     83      !!      ========  (excluding the vertical flux proportional to dk[t] ) 
     84      !!      vertical fluxes associated with the rotated lateral mixing: 
     85      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ] 
     86      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  } 
     87      !!      take the horizontal divergence of the fluxes: 
     88      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ] 
     89      !!      Add this trend to the general trend (ta,sa): 
     90      !!         pta = pta + difft 
     91      !! 
     92      !! ** Action :   Update pta arrays with the before rotated diffusion 
     93      !!---------------------------------------------------------------------- 
     94      USE oce,   zftu => ua   ! use ua as workspace 
     95      USE oce,   zftv => va   ! use va as workspace 
     96      !! 
     97      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index 
     98      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
     99      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
     100      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
     101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
     102      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
     103      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef 
     104      !! 
     105      INTEGER  ::  ji, jj, jk,jn   ! dummy loop indices 
     106      INTEGER  ::  ip,jp,kp        ! dummy loop indices 
     107      INTEGER  ::  ierr            ! temporary integer 
     108      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
     109      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
     110      REAL(wp) ::  zcoef0, zbtr                  !   -      - 
     111      REAL(wp), DIMENSION(jpi,jpj,0:1) ::   zdkt               ! 2D+1 workspace 
     112      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw   ! 3D workspace 
    114113      ! 
    115     REAL(wp) ::   zslope_skew,zslope_iso,zslope2, zbu, zbv 
    116     REAL(wp) ::   ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt 
    117     REAL(wp) ::   ah,zah_slp,zaei_slp 
     114      REAL(wp) ::   zslope_skew,zslope_iso,zslope2, zbu, zbv 
     115      REAL(wp) ::   ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt 
     116      REAL(wp) ::   ah,zah_slp,zaei_slp 
    118117#if defined key_diaar5 
    119     REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                ! 2D workspace 
    120     REAL(wp)                         ::   zztmp              ! local scalar 
     118      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                ! 2D workspace 
     119      REAL(wp)                         ::   zztmp              ! local scalar 
    121120#endif 
    122121      !!---------------------------------------------------------------------- 
    123122 
    124     IF( kt == nit000 )  THEN 
    125        IF(lwp) WRITE(numout,*) 
    126        IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 
    127        IF(lwp) WRITE(numout,*) '                   WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 
    128        IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    129        ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 
    130        IF( ierr > 0 ) THEN 
    131           CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' )   ;   RETURN 
    132        ENDIF 
    133        IF( ln_traldf_gdia ) THEN 
    134           ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 
    135           IF( ierr > 0 ) THEN 
    136              CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' )   ;   RETURN 
    137           ENDIF 
    138        ENDIF 
    139     ENDIF 
    140  
    141       !                                                
     123      IF( kt == nit000 )  THEN 
     124         IF(lwp) WRITE(numout,*) 
     125         IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 
     126         IF(lwp) WRITE(numout,*) '                   WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL' 
     127         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     128         ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr ) 
     129         IF( ierr > 0 ) THEN 
     130            CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' )   ;   RETURN 
     131         ENDIF 
     132         IF( ln_traldf_gdia ) THEN 
     133            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr ) 
     134            IF( ierr > 0 ) THEN 
     135               CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' )   ;   RETURN 
     136            ENDIF 
     137         ENDIF 
     138      ENDIF 
     139 
    142140      !!---------------------------------------------------------------------- 
    143141      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv  
    … …  
    146144!!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 
    147145 
    148     ah_wslp2(:,:,:) = 0.0 
    149     IF( ln_traldf_gdia ) THEN 
    150        psix_eiv(:,:,:) = 0.0 
    151        psiy_eiv(:,:,:) = 0.0 
    152     ENDIF 
    153  
    154     DO ip=0,1 
    155        DO kp=0,1 
    156           DO jk=1,jpkm1 
    157              DO jj=1,jpjm1 
    158                 DO ji=1,fs_jpim1 
    159                    ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 
    160                    zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 
    161                    ah = fsahtu(ji,jj,jk)                                  !  fsaht(ji+ip,jj,jk) 
    162                    zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    163                    zslope2=(zslope_skew & 
    164                         & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur & 
     146      ah_wslp2(:,:,:) = 0.0 
     147      IF( ln_traldf_gdia ) THEN 
     148         psix_eiv(:,:,:) = 0.0 
     149         psiy_eiv(:,:,:) = 0.0 
     150      ENDIF 
     151 
     152      DO ip=0,1 
     153         DO kp=0,1 
     154            DO jk=1,jpkm1 
     155               DO jj=1,jpjm1 
     156                  DO ji=1,fs_jpim1 
     157                     ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp) 
     158                     zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 
     159                     ah = fsahtu(ji,jj,jk)                                  !  fsaht(ji+ip,jj,jk) 
     160                     zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
     161                     zslope2 = ( zslope_skew & 
     162                        &      - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur )**2 
     163                     ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 
     164                        & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 
     165                     IF( ln_traldf_gdia ) THEN 
     166                        zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 
     167                        psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 
     168                     ENDIF 
     169                  END DO 
     170               END DO 
     171            END DO 
     172         END DO 
     173      END DO 
     174      ! 
     175      DO jp=0,1 
     176         DO kp=0,1 
     177            DO jk=1,jpkm1 
     178               DO jj=1,jpjm1 
     179                  DO ji=1,fs_jpim1 
     180                     ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 
     181                     zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 
     182                     ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 
     183                     zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
     184                     zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 
    165185                        & )**2 
    166                    ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + & 
    167                         & ah*((zbu*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*zslope2 
    168                    IF( ln_traldf_gdia ) THEN 
    169                       zaei_slp = fsaeiw(ji+ip,jj,jk)*zslope_skew!fsaeit(ji+ip,jj,jk)*zslope_skew 
    170                       psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 
    171                    ENDIF 
    172                 END DO 
    173              END DO 
    174           END DO 
    175        END DO 
    176     END DO 
    177  
    178     DO jp=0,1 
    179        DO kp=0,1 
    180           DO jk=1,jpkm1 
    181              DO jj=1,jpjm1 
    182                 DO ji=1,fs_jpim1 
    183                    ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp) 
    184                    zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 
    185                    ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk) 
    186                    zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    187                    zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr & 
    188                         & )**2 
    189                    ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 
     186                     ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + & 
    190187                        & ah*((zbv*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*zslope2 
    191                    IF( ln_traldf_gdia ) THEN 
    192                       zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew     !fsaeit(ji,jj+jp,jk)*zslope_skew 
    193                       psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 
    194                    ENDIF 
    195                 END DO 
    196              END DO 
    197           END DO 
    198        END DO 
     188                     IF( ln_traldf_gdia ) THEN 
     189                        zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew     !fsaeit(ji,jj+jp,jk)*zslope_skew 
     190                        psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zaei_slp 
     191                     ENDIF 
     192                  END DO 
     193               END DO 
     194            END DO 
     195         END DO 
    199196      END DO 
    200  
    201197      ! 
    202198      !                                                          ! =========== 
    … …  
    224220               DO ji = 1, jpim1 
    225221# endif 
    226                   iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1                ! last level 
    227                   ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
    228                   zdit(ji,jj,iku) = pgu(ji,jj,jn)           
    229                   zdjt(ji,jj,ikv) = pgv(ji,jj,jn)       
     222                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     223                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)       
    230224               END DO 
    231225            END DO 
    … …  
    294288            END DO 
    295289 
    296           !                        !==  divergence and add to the general trend  ==! 
    297           DO jj = 2 , jpjm1 
    298              DO ji = fs_2, fs_jpim1   ! vector opt. 
    299                 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    300                 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (   zftu(ji-1,jj,jk) - zftu(ji,jj,jk)   & 
    301                   &                                            + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   ) 
    302              END DO 
    303           END DO 
    304           ! 
    305        END DO 
    306        ! 
    307        DO jk = 1, jpkm1            !== Divergence of vertical fluxes added to the general tracer trend 
    308           DO jj = 2, jpjm1 
    309              DO ji = fs_2, fs_jpim1   ! vector opt. 
    310                 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
    311                   &                                 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    312              END DO 
    313           END DO 
    314        END DO 
    315        ! 
    316        !                            ! "Poleward" diffusive heat or salt transports (T-S case only) 
    317        IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
    318          IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) )        ! 3.3  names 
    319          IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) ) 
    320        ENDIF 
     290            !                        !==  divergence and add to the general trend  ==! 
     291            DO jj = 2 , jpjm1 
     292               DO ji = fs_2, fs_jpim1   ! vector opt. 
     293                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     294                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (   zftu(ji-1,jj,jk) - zftu(ji,jj,jk)   & 
     295                     &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   ) 
     296               END DO 
     297            END DO 
     298            ! 
     299         END DO 
     300         ! 
     301         DO jk = 1, jpkm1            !== Divergence of vertical fluxes added to the general tracer trend 
     302            DO jj = 2, jpjm1 
     303               DO ji = fs_2, fs_jpim1   ! vector opt. 
     304                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     305                     &                                / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     306               END DO 
     307            END DO 
     308         END DO 
     309         ! 
     310         !                            ! "Poleward" diffusive heat or salt transports (T-S case only) 
     311         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
     312            IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) )        ! 3.3  names 
     313            IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) ) 
     314         ENDIF 
    321315 
    322316#if defined key_diaar5 
    323        IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    324           zztmp = 0.5 * rau0 * rcp  
    325           z2d(:,:) = 0.e0  
    326           DO jk = 1, jpkm1 
    327              DO jj = 2, jpjm1 
    328                 DO ji = fs_2, fs_jpim1   ! vector opt. 
    329                    z2d(ji,jj) = z2d(ji,jj) + zztmp * zftu(ji,jj,jk)       & 
    330                         &       * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk)  
    331                 END DO 
    332              END DO 
    333           END DO 
    334           CALL lbc_lnk( z2d, 'U', -1. ) 
    335           CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    336           z2d(:,:) = 0.e0  
    337           DO jk = 1, jpkm1 
    338              DO jj = 2, jpjm1 
    339                 DO ji = fs_2, fs_jpim1   ! vector opt. 
    340                    z2d(ji,jj) = z2d(ji,jj) + zztmp * zftv(ji,jj,jk)       & 
    341                         &       * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk)  
    342                 END DO 
    343              END DO 
    344           END DO 
    345           CALL lbc_lnk( z2d, 'V', -1. ) 
    346           CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
    347        END IF 
     317         IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     318            z2d(:,:) = 0._wp  
     319            zztmp = rau0 * rcp  
     320            DO jk = 1, jpkm1 
     321               DO jj = 2, jpjm1 
     322                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     323                     z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     324                  END DO 
     325               END DO 
     326            END DO 
     327            z2d(:,:) = zztmp * z2d(:,:) 
     328            CALL lbc_lnk( z2d, 'U', -1. ) 
     329            CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     330            DO jk = 1, jpkm1 
     331               DO jj = 2, jpjm1 
     332                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     333                     z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     334                  END DO 
     335               END DO 
     336            END DO 
     337            z2d(:,:) = zztmp * z2d(:,:) 
     338            CALL lbc_lnk( z2d, 'V', -1. ) 
     339            CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
     340         END IF 
    348341#endif 
    349        ! 
    350     END DO 
    351     ! 
     342         ! 
     343      END DO 
     344      ! 
    352345  END SUBROUTINE tra_ldf_iso_grif 
    353346 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r2399 r2450  
    104104                  DO ji = 1, fs_jpim1   ! vector opt. 
    105105                     ! last level 
    106                      iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1 
    107                      ikv = MIN ( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
     106                     iku = mbku(ji,jj) 
     107                     ikv = mbkv(ji,jj) 
    108108                     IF( iku == jk ) THEN 
    109109                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku) 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r2287 r2450  
    3333   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3434   !! $Id$  
    35    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    36    !!---------------------------------------------------------------------- 
    37  
     35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     36   !!---------------------------------------------------------------------- 
    3837CONTAINS 
    3938 
    … …  
    131130                  IF( zwy(ji,1) /= 0.e0 ) THEN 
    132131                     ! 
    133                      ikbot = mbathy(ji,jj)      ! ikbot: ocean bottom level 
     132                     ikbot = mbkt(ji,jj)        ! ikbot: ocean bottom T-level 
    134133                     ! 
    135134                     DO jiter = 1, jpk          ! vertical iteration 
    … …  
    140139220                     CONTINUE 
    141140                        ik = ik + 1 
    142                         IF( ik >= ikbot-1 ) GO TO 200 
     141                        IF( ik >= ikbot ) GO TO 200 
    143142                        zwx(ji,ik) = zrhop(ji,jj,ik) - zrhop(ji,jj,ik+1) 
    144143                        IF( zwx(ji,ik) <= 0.e0 ) GO TO 220 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r2287 r2450  
    11MODULE zpshde 
    2    !!============================================================================== 
     2   !!====================================================================== 
    33   !!                       ***  MODULE zpshde   *** 
    4    !! z-coordinate - partial step : Horizontal Derivative 
    5    !!============================================================================== 
     4   !! z-coordinate + partial step : Horizontal Derivative at ocean bottom level 
     5   !!====================================================================== 
    66   !! History :  OPA  !  2002-04  (A. Bozec)  Original code 
    7    !!            8.5  !  2002-08  (G. Madec E. Durand)  Optimization and Free form 
    8    !!   NEMO     1.0  !  2004-03  (C. Ethe)  adapted for passive tracers 
     7   !!   NEMO     1.0  !  2002-08  (G. Madec E. Durand)  Optimization and Free form 
     8   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers 
    99   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA  
    10    !!============================================================================== 
     10   !!====================================================================== 
    1111    
    1212   !!---------------------------------------------------------------------- 
    … …  
    1414   !!                   ocean level (Z-coord. with Partial Steps) 
    1515   !!---------------------------------------------------------------------- 
    16    USE dom_oce         ! ocean space domain variables 
    17    USE oce             ! ocean dynamics and tracers variables 
     16   USE oce             ! ocean: dynamics and tracers variables 
     17   USE dom_oce         ! domain: ocean variables 
    1818   USE phycst          ! physical constants 
     19   USE eosbn2          ! ocean equation of state 
    1920   USE in_out_manager  ! I/O manager 
    20    USE eosbn2          ! ocean equation of state 
    2121   USE lbclnk          ! lateral boundary conditions (or mpp link) 
    2222 
    … …  
    2424   PRIVATE 
    2525 
    26    PUBLIC   zps_hde        ! routine called by step.F90 
    27    PUBLIC   zps_hde_init   ! routine called by opa.F90 
    28  
    29    INTEGER, DIMENSION(jpi,jpj) ::   mbatu, mbatv   ! bottom ocean level index at U- and V-points 
     26   PUBLIC   zps_hde    ! routine called by step.F90 
    3027 
    3128   !! * Substitutions 
    … …  
    3532   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    3633   !! $Id$ 
    37    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     34   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3835   !!---------------------------------------------------------------------- 
    3936CONTAINS 
    … …  
    4441      !!                     ***  ROUTINE zps_hde  *** 
    4542      !!                     
    46       !! ** Purpose :   Compute the horizontal derivative of T, S and rd 
     43      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
    4744      !!      at u- and v-points with a linear interpolation for z-coordinate 
    4845      !!      with partial steps. 
    … …  
    7875      !!      formulation of the equation of state (eos). 
    7976      !!      Gradient formulation for rho : 
    80       !!          di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~ 
    81       !! 
    82       !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at U/V-points 
    83       !!              - pgru, pgrv: horizontal gradient of rd if present at U/V-points  
    84       !!                and rd at V-points  
     77      !!          di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 
     78      !! 
     79      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
     80      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points  
    8581      !!---------------------------------------------------------------------- 
    8682      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    … …  
    9288      !! 
    9389      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    94       INTEGER  ::   iku, ikv        ! partial step level at u- and v-points 
     90      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    9591      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj     ! interpolated value of tracer 
    9692      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj     ! interpolated value of rd 
    … …  
    9995      !!---------------------------------------------------------------------- 
    10096 
    101  
    102       ! Interpolation of tracers at the last ocean level 
    103       DO jn = 1, kjpt 
     97      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    10498         ! 
    10599# if defined key_vectopt_loop 
    … …  
    110104            DO ji = 1, jpim1 
    111105# endif 
    112                ! last level 
    113                iku = mbatu(ji,jj) 
    114                ikv = mbatv(ji,jj) 
    115                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    116                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    117  
     106               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku , 1 )        ! last and before last ocean level at u- & v-points 
     107               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv , 1 )        ! if level first is a p-step, ik.m1=1 
     108               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     109               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     110               ! 
    118111               ! i- direction 
    119                IF( ze3wu >= 0. ) THEN      ! case 1 
     112               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    120113                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    121114                  ! interpolated values of tracers 
    122                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku-1,jn) - pta(ji+1,jj,iku,jn) ) 
     115                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
    123116                  ! gradient of  tracers 
    124117                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    125                ELSE                        ! case 2 
     118               ELSE                           ! case 2 
    126119                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    127120                  ! interpolated values of tracers 
    128                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku-1,jn) - pta(ji,jj,iku,jn) ) 
     121                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    129122                  ! gradient of tracers 
    130123                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    131124               ENDIF 
    132  
     125               ! 
    133126               ! j- direction 
    134                IF( ze3wv >= 0. ) THEN      ! case 1 
     127               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    135128                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    136129                  ! interpolated values of tracers 
    137                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv-1,jn) - pta(ji,jj+1,ikv,jn) ) 
     130                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    138131                  ! gradient of tracers 
    139132                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    140                ELSE                        ! case 2 
     133               ELSE                           ! case 2 
    141134                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    142135                  ! interpolated values of tracers 
    143                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv-1,jn) - pta(ji,jj,ikv,jn) ) 
     136                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    144137                  ! gradient of tracers 
    145138                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    … …  
    162155            DO ji = 1, jpim1 
    163156# endif 
    164                iku = mbatu(ji,jj) 
    165                ikv = mbatv(ji,jj) 
     157               iku = mbku(ji,jj) 
     158               ikv = mbkv(ji,jj) 
    166159               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    167160               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    168                IF( ze3wu >= 0. ) THEN    
    169                   zhi(ji,jj) = fsdept(ji  ,jj,iku) 
    170                ELSE                    
    171                   zhi(ji,jj) = fsdept(ji+1,jj,iku) 
    172                ENDIF 
    173                IF( ze3wv >= 0. ) THEN   
    174                   zhj(ji,jj) = fsdept(ji,jj  ,ikv) 
    175                ELSE                    
    176                   zhj(ji,jj) = fsdept(ji,jj+1,ikv) 
     161               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
     162               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
     163               ENDIF 
     164               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
     165               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
    177166               ENDIF 
    178167# if ! defined key_vectopt_loop 
    … …  
    193182            DO ji = 1, jpim1 
    194183# endif 
    195                iku = mbatu(ji,jj) 
    196                ikv = mbatv(ji,jj) 
     184               iku = mbku(ji,jj) 
     185               ikv = mbkv(ji,jj) 
    197186               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    198187               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    199                IF( ze3wu >= 0. ) THEN    ! i-direction: case 1 
    200                   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) ) 
    201                ELSE                      ! i-direction: case 2 
    202                   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) ) 
    203                ENDIF 
    204                IF( ze3wv >= 0. ) THEN    ! j-direction: case 1 
    205                   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) ) 
    206                ELSE                      ! j-direction: case 2 
    207                   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) ) 
     188               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
     189               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
     190               ENDIF 
     191               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
     192               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    208193               ENDIF 
    209194# if ! defined key_vectopt_loop 
    … …  
    217202   END SUBROUTINE zps_hde 
    218203 
    219  
    220    SUBROUTINE zps_hde_init 
    221       !!---------------------------------------------------------------------- 
    222       !!                     ***  ROUTINE zps_hde_init  *** 
    223       !! 
    224       !! ** Purpose : Computation of bottom ocean level index at U- and V-points  
    225       !!                     
    226       !!---------------------------------------------------------------------- 
    227       INTEGER ::   ji, jj   ! Dummy loop indices 
    228       REAL(wp), DIMENSION(jpi,jpj) ::   zti, ztj     ! 2D workspace  
    229       !!---------------------------------------------------------------------- 
    230       ! 
    231       mbatu(:,:) = 0 
    232       mbatv(:,:) = 0 
    233       DO jj = 1, jpjm1 
    234          DO ji = 1, fs_jpim1   ! vector opt. 
    235             mbatu(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1, 2 ) 
    236             mbatv(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1, 2 ) 
    237          END DO 
    238       END DO 
    239       zti(:,:) = FLOAT( mbatu(:,:) ) 
    240       ztj(:,:) = FLOAT( mbatv(:,:) ) 
    241       ! lateral boundary conditions: T-point, sign unchanged 
    242       CALL lbc_lnk( zti , 'U', 1. )   ;   CALL lbc_lnk( ztj , 'V', 1. ) 
    243       mbatu(:,:) = MAX( INT( zti(:,:) ), 2 ) 
    244       mbatv(:,:) = MAX( INT( ztj(:,:) ), 2 ) 
    245       ! 
    246    END SUBROUTINE zps_hde_init 
    247204   !!====================================================================== 
    248205END MODULE zpshde 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdvor.F90

    r2364 r2450  
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_trdvor'   : momentum trend diagnostics 
    12    !!---------------------------------------------------------------------- 
    1312   !!---------------------------------------------------------------------- 
    1413   !!   trd_vor      : momentum trends averaged over the depth 
    … …  
    3938   PUBLIC   trd_vor_init   ! routine called by opa.F90 
    4039 
    41    INTEGER ::                & 
    42       nh_t, nmoydpvor  ,     & 
    43       nidvor, nhoridvor,     & 
    44       ndexvor1(jpi*jpj),     & 
    45       ndimvor1, icount,      & 
    46       idebug                    ! (0/1) set it to 1 in case of problem to have more print 
    47  
    48    REAL(wp), DIMENSION(jpi,jpj) ::  & 
    49      vor_avr    ,     &  ! average 
    50      vor_avrb   ,     &  ! before vorticity (kt-1) 
    51      vor_avrbb  ,     &  ! vorticity at begining of the nwrite-1 timestep averaging period 
    52      vor_avrbn  ,     &  ! after vorticity at time step after the 
    53      rotot      ,     &  ! begining of the NWRITE-1 timesteps 
    54      vor_avrtot ,     & 
    55      vor_avrres 
    56  
    57    REAL(wp), DIMENSION(jpi,jpj,jpltot_vor)::   vortrd  !: curl of trends 
     40   INTEGER ::   nh_t, nmoydpvor, nidvor, nhoridvor, ndexvor1(jpi*jpj), ndimvor1, icount   ! needs for IOIPSL output 
     41   INTEGER ::   ndebug     ! (0/1) set it to 1 in case of problem to have more print 
     42 
     43   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avr      ! average 
     44   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrb     ! before vorticity (kt-1) 
     45   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrbb    ! vorticity at begining of the nwrite-1 timestep averaging period 
     46   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrbn    ! after vorticity at time step after the 
     47   REAL(wp), DIMENSION(jpi,jpj) ::   rotot        ! begining of the NWRITE-1 timesteps 
     48   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrtot   ! 
     49   REAL(wp), DIMENSION(jpi,jpj) ::   vor_avrres   ! 
     50 
     51   REAL(wp), DIMENSION(jpi,jpj,jpltot_vor) ::   vortrd  ! curl of trends 
    5852          
    5953   CHARACTER(len=12) ::   cvort 
    … …  
    6660   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    6761   !! $Id$  
    68    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    69    !!---------------------------------------------------------------------- 
    70    
     62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     63   !!---------------------------------------------------------------------- 
    7164CONTAINS 
    7265 
    … …  
    7669      !! 
    7770      !! ** Purpose :   computation of vertically integrated vorticity budgets 
    78       !!      from ocean surface down to control surface (NetCDF output) 
    79       !! 
    80       !! ** Method/usage : 
    81       !!      integration done over nwrite-1 time steps 
    82       !! 
    83       !! 
    84       !! ** Action : 
    85       !!            /comvor/   : 
    86       !!                         vor_avr          average 
    87       !!                         vor_avrb         vorticity at kt-1 
    88       !!                         vor_avrbb        vorticity at begining of the NWRITE-1 
    89       !!                                          time steps averaging period 
    90       !!                         vor_avrbn         vorticity at time step after the 
    91       !!                                          begining of the NWRITE-1 time 
    92       !!                                          steps averaging period 
    93       !! 
    94       !!                 trends : 
    95       !! 
     71      !!              from ocean surface down to control surface (NetCDF output) 
     72      !! 
     73      !! ** Method/usage :   integration done over nwrite-1 time steps 
     74      !! 
     75      !! ** Action :   trends : 
    9676      !!                  vortrd (,, 1) = Pressure Gradient Trend 
    9777      !!                  vortrd (,, 2) = KE Gradient Trend 
    … …  
    10484      !!                  vortrd (,, 9) = Beta V 
    10585      !!                  vortrd (,,10) = forcing term 
    106       !!      vortrd (,,11) = bottom friction term 
     86      !!                  vortrd (,,11) = bottom friction term 
    10787      !!                  rotot(,) : total cumulative trends over nwrite-1 time steps 
    10888      !!                  vor_avrtot(,) : first membre of vrticity equation 
    … …  
    11090      !! 
    11191      !!      trends output in netCDF format using ioipsl 
    112       !! 
    113       !!---------------------------------------------------------------------- 
    114       INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index 
    115       REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   & 
    116          putrdvor,                         &  ! u vorticity trend  
    117          pvtrdvor                             ! v vorticity trend 
    118       !! 
    119       INTEGER ::   ji, jj 
    120       INTEGER ::   ikbu, ikbum1, ikbv, ikbvm1 
     92      !!---------------------------------------------------------------------- 
     93      INTEGER                     , INTENT(in   ) ::   ktrd       ! ocean trend index 
     94      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   putrdvor   ! u vorticity trend  
     95      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pvtrdvor   ! v vorticity trend 
     96      !! 
     97      INTEGER ::   ji, jj       ! dummy loop indices 
     98      INTEGER ::   ikbu, ikbv   ! local integers 
    12199      REAL(wp), DIMENSION(jpi,jpj) ::   zudpvor, zvdpvor   ! total cmulative trends 
    122100      !!---------------------------------------------------------------------- 
    123101 
    124102      ! Initialization 
    125       zudpvor(:,:) = 0.e0 
    126       zvdpvor(:,:) = 0.e0 
    127  
    128       CALL lbc_lnk( putrdvor,  'U' , -1. ) 
     103      zudpvor(:,:) = 0._wp 
     104      zvdpvor(:,:) = 0._wp 
     105      ! 
     106      CALL lbc_lnk( putrdvor,  'U' , -1. )         ! lateral boundary condition on input momentum trends 
    129107      CALL lbc_lnk( pvtrdvor,  'V' , -1. ) 
    130108 
    … …  
    134112 
    135113      SELECT CASE (ktrd)  
    136  
     114      ! 
    137115      CASE (jpvor_bfr)        ! bottom friction 
    138  
    139116         DO jj = 2, jpjm1 
    140117            DO ji = fs_2, fs_jpim1  
    141                ikbu   = min( mbathy(ji+1,jj), mbathy(ji,jj) ) 
    142                ikbum1 = max( ikbu-1, 1 ) 
    143                ikbv   = min( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    144                ikbvm1 = max( ikbv-1, 1 ) 
    145              
    146                zudpvor(ji,jj) = putrdvor(ji,jj) * fse3u(ji,jj,ikbum1) * e1u(ji,jj) * umask(ji,jj,ikbum1) 
    147                zvdpvor(ji,jj) = pvtrdvor(ji,jj) * fse3v(ji,jj,ikbvm1) * e2v(ji,jj) * vmask(ji,jj,ikbvm1) 
     118               ikbu = mbkv(ji,jj) 
     119               ikbv = mbkv(ji,jj)             
     120               zudpvor(ji,jj) = putrdvor(ji,jj) * fse3u(ji,jj,ikbu) * e1u(ji,jj) * umask(ji,jj,ikbu) 
     121               zvdpvor(ji,jj) = pvtrdvor(ji,jj) * fse3v(ji,jj,ikbv) * e2v(ji,jj) * vmask(ji,jj,ikbv) 
    148122            END DO 
    149123         END DO 
    150  
     124         ! 
    151125      CASE (jpvor_swf)        ! wind stress 
    152  
    153126         zudpvor(:,:) = putrdvor(:,:) * fse3u(:,:,1) * e1u(:,:) * umask(:,:,1) 
    154127         zvdpvor(:,:) = pvtrdvor(:,:) * fse3v(:,:,1) * e2v(:,:) * vmask(:,:,1) 
    155  
     128         ! 
    156129      END SELECT 
    157130 
    … …  
    163136      DO ji=1,jpim1 
    164137         DO jj=1,jpjm1 
    165             vortrd(ji,jj,ktrd) = ( zvdpvor(ji+1,jj) - zvdpvor(ji,jj)        & 
    166                  &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) & 
    167                  &               / ( e1f(ji,jj) * e2f(ji,jj) ) 
     138            vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)       & 
     139                 &                - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) )   ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    168140         END DO 
    169141      END DO 
    170  
    171       ! Surface mask 
    172       vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1) 
    173  
    174       IF( idebug /= 0 ) THEN 
     142      vortrd(:,:,ktrd) = vortrd(:,:,ktrd) * fmask(:,:,1)      ! Surface mask 
     143 
     144      IF( ndebug /= 0 ) THEN 
    175145         IF(lwp) WRITE(numout,*) ' debuging trd_vor_zint: I done' 
    176146         CALL FLUSH(numout) 
    … …  
    185155      !! 
    186156      !! ** Purpose :   computation of vertically integrated vorticity budgets 
    187       !!      from ocean surface down to control surface (NetCDF output) 
    188       !! 
    189       !! ** Method/usage : 
    190       !!      integration done over nwrite-1 time steps 
    191       !! 
    192       !! 
    193       !! ** Action : 
    194       !!            /comvor/   : 
    195       !!                         vor_avr          average 
    196       !!                         vor_avrb         vorticity at kt-1 
    197       !!                         vor_avrbb        vorticity at begining of the NWRITE-1 
    198       !!                                          time steps averaging period 
    199       !!                         vor_avrbn         vorticity at time step after the 
    200       !!                                          begining of the NWRITE-1 time 
    201       !!                                          steps averaging period 
    202       !! 
    203       !!                 trends : 
    204       !! 
     157      !!              from ocean surface down to control surface (NetCDF output) 
     158      !! 
     159      !! ** Method/usage :   integration done over nwrite-1 time steps 
     160      !! 
     161      !! ** Action :     trends : 
    205162      !!                  vortrd (,,1) = Pressure Gradient Trend 
    206163      !!                  vortrd (,,2) = KE Gradient Trend 
    … …  
    219176      !! 
    220177      !!      trends output in netCDF format using ioipsl 
    221       !! 
    222       !!---------------------------------------------------------------------- 
    223       INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index 
    224       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   & 
    225          putrdvor,                         &  ! u vorticity trend  
    226          pvtrdvor                             ! v vorticity trend 
     178      !!---------------------------------------------------------------------- 
     179      INTEGER                         , INTENT(in   ) ::   ktrd       ! ocean trend index 
     180      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   putrdvor   ! u vorticity trend  
     181      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvtrdvor   ! v vorticity trend 
    227182      !! 
    228183      INTEGER ::   ji, jj, jk 
    229       REAL(wp), DIMENSION(jpi,jpj) ::   & 
    230          zubet,                         &  ! u Beta.V case 
    231          zvbet,                         &  ! v Beta.V case 
    232          zudpvor,                       &  ! total cmulative trends 
    233          zvdpvor                           !   "      "        " 
     184      REAL(wp), DIMENSION(jpi,jpj) ::   zubet  , zvbet     ! Beta.V  
     185      REAL(wp), DIMENSION(jpi,jpj) ::   zudpvor, zvdpvor   ! total cmulative trends 
    234186      !!---------------------------------------------------------------------- 
    235187      
    236188      ! Initialization 
    237       zubet(:,:) = 0.e0 
    238       zvbet(:,:) = 0.e0 
    239       zudpvor(:,:) = 0.e0 
    240       zvdpvor(:,:) = 0.e0 
     189      zubet  (:,:) = 0._wp 
     190      zvbet  (:,:) = 0._wp 
     191      zudpvor(:,:) = 0._wp 
     192      zvdpvor(:,:) = 0._wp 
     193      ! 
     194      CALL lbc_lnk( putrdvor, 'U' , -1. )         ! lateral boundary condition on input momentum trends 
     195      CALL lbc_lnk( pvtrdvor, 'V' , -1. ) 
    241196 
    242197      !  ===================================== 
    243198      !  I vertical integration of 3D trends 
    244199      !  ===================================== 
    245  
    246       CALL lbc_lnk( putrdvor, 'U' , -1. ) 
    247       CALL lbc_lnk( pvtrdvor, 'V' , -1. ) 
    248  
    249200      ! putrdvor and pvtrdvor terms 
    250201      DO jk = 1,jpk 
    … …  
    267218      DO ji=1,jpim1 
    268219         DO jj=1,jpjm1 
    269             vortrd(ji,jj,ktrd) = (  zvdpvor(ji+1,jj) - zvdpvor(ji,jj) -   & 
    270                  &                ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) & 
    271                  &               / ( e1f(ji,jj) * e2f(ji,jj) ) 
     220            vortrd(ji,jj,ktrd) = (    zvdpvor(ji+1,jj) - zvdpvor(ji,jj)     & 
     221               &                  - ( zudpvor(ji,jj+1) - zudpvor(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    272222         END DO 
    273223      END DO 
    … …  
    281231         DO ji=1,jpim1 
    282232            DO jj=1,jpjm1 
    283                vortrd(ji,jj,jpvor_bev) = (  zvbet(ji+1,jj) - zvbet(ji,jj) -   & 
    284                     &                    ( zubet(ji,jj+1) - zubet(ji,jj) ) ) & 
    285                     &                   / ( e1f(ji,jj) * e2f(ji,jj) ) 
     233               vortrd(ji,jj,jpvor_bev) = (    zvbet(ji+1,jj) - zvbet(ji,jj)     & 
     234                  &                       - ( zubet(ji,jj+1) - zubet(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) 
    286235            END DO 
    287236         END DO 
    … …  
    294243      ENDIF 
    295244    
    296       IF( idebug /= 0 ) THEN 
     245      IF( ndebug /= 0 ) THEN 
    297246         IF(lwp) WRITE(numout,*) ' debuging trd_vor_zint: I done' 
    298247         CALL FLUSH(numout) 
    … …  
    309258      !!               and make outputs (NetCDF or DIMG format) 
    310259      !!---------------------------------------------------------------------- 
    311       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    312       !! 
    313       INTEGER  ::   ji, jj, jk, jl, it, itmod 
    314       REAL(wp) ::   zmean 
    315       REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn 
     260      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     261      !! 
     262      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     263      INTEGER  ::   it, itmod        ! local integers 
     264      REAL(wp) ::   zmean            ! local scalars 
     265      REAL(wp), DIMENSION(jpi,jpj) ::   zun, zvn   ! 2D workspace 
    316266      !!---------------------------------------------------------------------- 
    317267 
    … …  
    324274      ! --------------------------------------------------- 
    325275 
    326       IF( kt > nit000 ) THEN 
    327          vor_avrb(:,:) = vor_avr(:,:) 
    328       ENDIF 
    329  
    330        IF( idebug /= 0 ) THEN 
     276      IF( kt > nit000 )   vor_avrb(:,:) = vor_avr(:,:) 
     277 
     278      IF( ndebug /= 0 ) THEN 
    331279          WRITE(numout,*) ' debuging trd_vor: I.1 done ' 
    332280          CALL FLUSH(numout) 
    … …  
    336284      !  ---------------------------------- 
    337285 
    338       vor_avr(:,:) = 0. 
    339       zun(:,:)=0 
    340       zvn(:,:)=0 
    341       vor_avrtot(:,:)=0 
    342       vor_avrres(:,:)=0 
     286      vor_avr   (:,:) = 0._wp 
     287      zun       (:,:) = 0._wp 
     288      zvn       (:,:) = 0._wp 
     289      vor_avrtot(:,:) = 0._wp 
     290      vor_avrres(:,:) = 0._wp 
    343291       
    344292      ! Vertically averaged velocity 
    345293      DO jk = 1, jpk - 1 
    346          zun(:,:)=zun(:,:) + e1u(:,:)*un(:,:,jk)*fse3u(:,:,jk) 
    347          zvn(:,:)=zvn(:,:) + e2v(:,:)*vn(:,:,jk)*fse3v(:,:,jk) 
     294         zun(:,:) = zun(:,:) + e1u(:,:) * un(:,:,jk) * fse3u(:,:,jk) 
     295         zvn(:,:) = zvn(:,:) + e2v(:,:) * vn(:,:,jk) * fse3v(:,:,jk) 
    348296      END DO 
    349297  
    350       zun(:,:)=zun(:,:)*hur(:,:) 
    351       zvn(:,:)=zvn(:,:)*hvr(:,:) 
     298      zun(:,:) = zun(:,:) * hur(:,:) 
     299      zvn(:,:) = zvn(:,:) * hvr(:,:) 
    352300 
    353301      ! Curl 
    354302      DO ji=1,jpim1 
    355303         DO jj=1,jpjm1 
    356             vor_avr(ji,jj) = ((zvn(ji+1,jj)-zvn(ji,jj))-   & 
    357                               (zun(ji,jj+1)-zun(ji,jj)))   & 
    358                              /( e1f(ji,jj) * e2f(ji,jj) ) 
    359             vor_avr(ji,jj) = vor_avr(ji,jj)*fmask(ji,jj,1) 
     304            vor_avr(ji,jj) = (  ( zvn(ji+1,jj) - zvn(ji,jj) )    & 
     305               &              - ( zun(ji,jj+1) - zun(ji,jj) ) ) / ( e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,1) 
    360306         END DO 
    361307      END DO 
    362308       
    363       IF(idebug /= 0) THEN 
     309      IF( ndebug /= 0 ) THEN 
    364310         WRITE(numout,*) ' debuging trd_vor: I.2 done' 
    365311         CALL FLUSH(numout) 
    … …  
    377323      ENDIF 
    378324 
    379       IF( idebug /= 0 ) THEN 
     325      IF( ndebug /= 0 ) THEN 
    380326         WRITE(numout,*) ' debuging trd_vor: I1.1 done' 
    381327         CALL FLUSH(numout) 
    … …  
    395341      ENDIF 
    396342 
    397       IF( idebug /= 0 ) THEN 
     343      IF( ndebug /= 0 ) THEN 
    398344         WRITE(numout,*) ' debuging trd_vor: II.2 done' 
    399345         CALL FLUSH(numout) 
    … …  
    405351 
    406352      ! define time axis 
    407       it = kt 
     353      it    = kt 
    408354      itmod = kt - nit000 + 1 
    409355 
    … …  
    412358         ! III.1 compute total trend 
    413359         ! ------------------------ 
    414          zmean = float(nmoydpvor) 
    415  
    416          vor_avrtot(:,:) = ( vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - & 
    417                              vor_avrbb(:,:) ) /  (zmean * 2. * rdt) 
    418  
    419          IF( idebug /= 0 ) THEN 
     360         zmean = 1._wp / (  REAL( nmoydpvor, wp ) * 2._wp * rdt  ) 
     361         vor_avrtot(:,:) = (  vor_avr(:,:) - vor_avrbn(:,:) + vor_avrb(:,:) - vor_avrbb(:,:) ) * zmean 
     362 
     363         IF( ndebug /= 0 ) THEN 
    420364             WRITE(numout,*) ' zmean = ',zmean 
    421365             WRITE(numout,*) ' debuging trd_vor: III.1 done' 
    … …  
    425369         ! III.2 compute residual 
    426370         ! --------------------- 
     371         zmean = 1._wp / REAL( nmoydpvor, wp ) 
    427372         vor_avrres(:,:) = vor_avrtot(:,:) - rotot(:,:) / zmean 
    428373 
    … …  
    431376         CALL lbc_lnk( vor_avrres, 'F', 1. ) 
    432377 
    433          IF( idebug /= 0 ) THEN 
     378         IF( ndebug /= 0 ) THEN 
    434379            WRITE(numout,*) ' debuging trd_vor: III.2 done' 
    435380            CALL FLUSH(numout) 
    … …  
    439384         ! ------------------------------ 
    440385         vor_avrbb(:,:) = vor_avrb(:,:) 
    441          vor_avrbn(:,:) = vor_avr(:,:) 
    442  
    443          IF( idebug /= 0 ) THEN 
     386         vor_avrbn(:,:) = vor_avr (:,:) 
     387 
     388         IF( ndebug /= 0 ) THEN 
    444389            WRITE(numout,*) ' debuging trd_vor: III.3 done' 
    445390            CALL FLUSH(numout) 
    446391         ENDIF 
    447  
    448          nmoydpvor=0 
    449  
     392         ! 
     393         nmoydpvor = 0 
     394         ! 
    450395      ENDIF 
    451396 
    … …  
    475420         CALL histwrite( nidvor,"sovorgap",it,vor_avrres    ,ndimvor1,ndexvor1) ! gap between 1st and 2 nd mbre 
    476421         ! 
    477          IF( idebug /= 0 ) THEN 
     422         IF( ndebug /= 0 ) THEN 
    478423            WRITE(numout,*) ' debuging trd_vor: III.4 done' 
    479424            CALL FLUSH(numout) 
    … …  
    508453 
    509454      ! Open specifier 
    510       idebug = 0      ! set it to 1 in case of problem to have more Print 
     455      ndebug = 0      ! set it to 1 in case of problem to have more Print 
    511456 
    512457      IF(lwp) THEN 
    … …  
    528473      vor_avrres(:,:)=0 
    529474 
    530       IF( idebug /= 0 ) THEN 
     475      IF( ndebug /= 0 ) THEN 
    531476         WRITE(numout,*) ' debuging trd_vor_init: I. done' 
    532477         CALL FLUSH(numout) 
    … …  
    600545      CALL histend( nidvor, snc4set ) 
    601546 
    602       IF( idebug /= 0 ) THEN 
     547      IF( ndebug /= 0 ) THEN 
    603548         WRITE(numout,*) ' debuging trd_vor_init: II. done' 
    604549         CALL FLUSH(numout) 
    … …  
    621566      REAL, DIMENSION(:,:), INTENT( inout ) ::   putrdvor, pvtrdvor 
    622567      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index 
    623       WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1) 
    624       WRITE(*,*) '  "      "     : You should not have seen this print! error?', pvtrdvor(1,1) 
    625       WRITE(*,*) '  "      "     : You should not have seen this print! error?', ktrd 
     568      WRITE(*,*) 'trd_vor_zint_2d: You should not have seen this print! error?', putrdvor(1,1), pvtrdvor(1,1), ktrd 
    626569   END SUBROUTINE trd_vor_zint_2d 
    627570   SUBROUTINE trd_vor_zint_3d( putrdvor, pvtrdvor, ktrd ) 
    628571      REAL, DIMENSION(:,:,:), INTENT( inout ) ::   putrdvor, pvtrdvor 
    629572      INTEGER, INTENT( in ) ::   ktrd         ! ocean trend index 
    630       WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1) 
    631       WRITE(*,*) '  "      "     : You should not have seen this print! error?', pvtrdvor(1,1,1) 
    632       WRITE(*,*) '  "      "     : You should not have seen this print! error?', ktrd 
     573      WRITE(*,*) 'trd_vor_zint_3d: You should not have seen this print! error?', putrdvor(1,1,1), pvtrdvor(1,1,1), ktrd 
    633574   END SUBROUTINE trd_vor_zint_3d 
    634575   SUBROUTINE trd_vor_init              ! Empty routine 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r2346 r2450  
    4040   REAL(wp) ::   rn_bfri2  = 1.0e-3_wp   ! bottom drag coefficient (non linear case) 
    4141   REAL(wp) ::   rn_bfeb2  = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2] 
    42    REAL(wp) ::   rn_bfrien = 30_wp       ! local factor to enhance coefficient bfri 
     42   REAL(wp) ::   rn_bfrien = 30._wp      ! local factor to enhance coefficient bfri 
    4343   LOGICAL  ::   ln_bfr2d  = .false.     ! logical switch for 2D enhancement 
    4444    
    … …  
    7474      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    7575      !! 
    76       INTEGER  ::   ji, jj         ! dummy loop indices 
    77       INTEGER  ::   ikbu, ikbum1   ! temporary integers 
    78       INTEGER  ::   ikbv, ikbvm1   !    -          - 
     76      INTEGER  ::   ji, jj       ! dummy loop indices 
     77      INTEGER  ::   ikbu, ikbv   ! local integers 
    7978      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars 
    8079      !!---------------------------------------------------------------------- 
    … …  
    9695            DO ji = 2, jpim1 
    9796# endif 
    98                ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    99                ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
    100                ikbum1 = MAX( ikbu-1, 1 ) 
    101                ikbvm1 = MAX( ikbv-1, 1 ) 
     97               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     98               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    10299               ! 
    103                zvu  = 0.25 * (  vn(ji,jj  ,ikbum1) + vn(ji+1,jj  ,ikbum1)     & 
    104                   &           + vn(ji,jj-1,ikbum1) + vn(ji+1,jj-1,ikbum1)  ) 
    105                zuv  = 0.25 * (  un(ji,jj  ,ikbvm1) + un(ji-1,jj  ,ikbvm1)     & 
    106                   &           + un(ji,jj+1,ikbvm1) + un(ji-1,jj+1,ikbvm1)  ) 
     100               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
     101                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
     102               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
     103                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
    107104               ! 
    108                zecu = SQRT(  un(ji,jj,ikbum1) * un(ji,jj,ikbum1) + zvu*zvu + rn_bfeb2  ) 
    109                zecv = SQRT(  vn(ji,jj,ikbvm1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2  ) 
     105               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  ) 
     106               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  ) 
    110107               ! 
    111                bfrua(ji,jj) = - 0.5 * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj  ) ) * zecu  
    112                bfrva(ji,jj) = - 0.5 * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji  ,jj+1) ) * zecv 
     108               bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj  ) ) * zecu  
     109               bfrva(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji  ,jj+1) ) * zecv 
    113110            END DO 
    114111         END DO 
    … …  
    126123      DO jj = 2, jpjm1 
    127124         DO ji = fs_2, fs_jpim1   ! vector opt. 
    128             !  
    129             ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    130             ikbum1 = MAX( ikbu-1, 1 ) 
    131             wbotu(ji,jj) = bfrua(ji,jj) * un(ji,jj,ikbum1) * umask(ji,jj,ikbum1) 
    132             ! 
    133             ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    134             ikbvm1 = MAX( ikbv-1, 1 ) 
    135             wbotv(ji,jj) = bfrva(ji,jj) * vn(ji,jj,ikbvm1) * vmask(ji,jj,ikbvm1) 
    136             ! 
     125            wbotu(ji,jj) = bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * umask(ji,jj,mbku(ji,jj)) 
     126            wbotv(ji,jj) = bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,mbkv(ji,jj)) 
    137127         END DO 
    138128      END DO 
    139       CALL lbc_lnk( wbotu(:,:), 'U', -1. ) 
    140       CALL lbc_lnk( wbotv(:,:), 'V', -1. ) 
     129      CALL lbc_lnk( wbotu(:,:), 'U', -1. )   ;   CALL lbc_lnk( wbotv(:,:), 'V', -1. )     ! lateral boundary condition 
    141130#endif 
    142131 
    … …  
    242231         DO ji = 2, jpim1 
    243232#  endif 
    244              ikbu = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    245              ikbv = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
     233             ikbu = mbku(ji,jj)       ! deepest ocean level at u- and v-points 
     234             ikbv = mbkv(ji,jj) 
    246235             zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 
    247236             zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 
    248237             IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN 
    249238                IF( ln_ctl ) THEN 
    250                    WRITE(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbu 
    251                    WRITE(numout,*) 'BFR ',ABS( bfrcoef2d(ji,jj) ), zfru 
     239                   WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu 
     240                   WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru 
    252241                ENDIF 
    253242                ictu = ictu + 1 
    … …  
    270259         CALL mpp_max( zmaxbfr ) 
    271260      ENDIF 
    272       IF( lwp .AND. ictu + ictv .GT. 0 ) THEN 
     261      IF( lwp .AND. ictu + ictv > 0 ) THEN 
    273262         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' 
    274263         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r2409 r2450  
    196196            DO jj = 2, jpjm1  
    197197               DO ji = fs_2, fs_jpim1   ! vector opt. 
    198                   zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbathy(ji,jj)) 
    199                   zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbathy(ji,jj)) ) 
     198                  zup   = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) 
     199                  zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) ) 
    200200                  zcoef = ( zup / MAX( zdown, rsmall ) ) 
    201201                  zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) 
    … …  
    360360!CDIR NOVERRCHK 
    361361            DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                ibot   = mbathy(ji,jj) 
    363                ibotm1 = ibot-1 
     362               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     363               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    364364               ! 
    365365               ! Bottom level Dirichlet condition: 
    … …  
    383383!CDIR NOVERRCHK 
    384384            DO ji = fs_2, fs_jpim1   ! vector opt. 
    385                ibot   = mbathy(ji,jj) 
    386                ibotm1 = ibot-1 
     385               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     386               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    387387               ! 
    388388               ! Bottom level Dirichlet condition: 
    … …  
    623623!CDIR NOVERRCHK 
    624624            DO ji = fs_2, fs_jpim1   ! vector opt. 
    625                ibot = mbathy(ji,jj) 
    626                ibotm1 = ibot-1 
     625               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     626               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    627627               zdep(ji,jj) = vkarmn * rn_hbro 
    628628               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    … …  
    646646!CDIR NOVERRCHK 
    647647            DO ji = fs_2, fs_jpim1   ! vector opt. 
    648                ibot   = mbathy(ji,jj) 
    649                ibotm1 = ibot-1 
     648               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
     649               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    650650               ! 
    651651               ! Bottom level Dirichlet condition: 
    … …  
    825825      DO jj = 2, jpjm1 
    826826         DO ji = fs_2, fs_jpim1   ! vector opt. 
    827             ibot   = mbathy(ji,jj) 
    828             avmv(ji,jj,ibot) = zcoef 
     827            avmv(ji,jj,mbkt(ji,jj)+1) = zcoef 
    829828         END DO 
    830829      END DO 
    … …  
    12161215              DO jj = 2, jpjm1 
    12171216                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1218                   ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    1219                   ikbum1 = MAX( ikbu-1, 1 ) 
    1220                   ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    1221                   ikbvm1 = MAX( ikbv-1, 1 ) 
     1217                  ikbu   = mbku(ji,jj) + 1    ! k   bottom level of uw-point 
     1218                  ikbum1 = mbku(ji,jj)        ! k-1 bottom level of u -point, but >=1 
     1219                  ikbv   = mbkv(ji,jj) + 1 
     1220                  ikbvm1 = mbkv(ji,jj) 
    12221221                  cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) 
    12231222                  cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) 
    1224                   wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1) 
    1225                   wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1) 
     1223                  wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1) * umask(ji,jj,1) 
     1224                  wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1) * vmask(ji,jj,1) 
    12261225                END DO 
    12271226              END DO 
    … …  
    12351234           ! Initialize bottom stresses 
    12361235           DO jj = 2, jpjm1 
    1237              DO ji = fs_2, fs_jpim1   ! vector opt. 
    1238                ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    1239                ikbum1 = MAX( ikbu-1, 1 ) 
    1240                ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    1241                ikbvm1 = MAX( ikbv-1, 1 ) 
    1242                cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) 
    1243                cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) 
    1244                wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1)*umask(ji,jj,1) 
    1245                wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1)*vmask(ji,jj,1) 
    1246              END DO 
     1236              DO ji = fs_2, fs_jpim1   ! vector opt. 
     1237                 ikbu   = mbku(ji,jj) + 1    ! k   bottom level of uw-point 
     1238                 ikbum1 = mbku(ji,jj)        ! k-1 bottom level of u -point, but >=1 
     1239                 ikbv   = mbkv(ji,jj) + 1 
     1240                 ikbvm1 = mbkv(ji,jj) 
     1241                 cbx = avmu(ji,jj,ikbu) / fse3uw(ji,jj,ikbu) 
     1242                 cby = avmv(ji,jj,ikbv) / fse3vw(ji,jj,ikbv) 
     1243                 wbotu(ji,jj) = -cbx * un(ji,jj,ikbum1) * umask(ji,jj,1) 
     1244                 wbotv(ji,jj) = -cby * vn(ji,jj,ikbvm1) * vmask(ji,jj,1) 
     1245              END DO 
    12471246           END DO 
    12481247        ENDIF 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90

    r2409 r2450  
    439439            zria(ji ) = 0. 
    440440            ! Maximum boundary layer depth 
    441             ikbot     = mbathy(ji,jj) - 1 ! ikbot is the last T point in the water 
     441            ikbot     = mbkt(ji,jj)    ! ikbot is the last T point in the water 
    442442            zhmax(ji) = fsdept(ji,jj,ikbot) - 0.001       
    443443            ! Compute Monin obukhov length scale at the surface and Ekman depth: 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r2287 r2450  
    6969 
    7070      ! w-level of the mixing and mixed layers 
    71       nmln(:,:) = mbathy(:,:)          ! Initialization to the number of w ocean point mbathy 
    72       imld(:,:) = mbathy(:,:) 
     71      nmln(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point 
     72      imld(:,:) = mbkt(:,:) + 1 
    7373      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
    7474         DO jj = 1, jpj 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r2386 r2450  
    225225!CDIR NOVERRCHK 
    226226!!       DO ji = fs_2, fs_jpim1   ! vector opt. 
    227 !!          ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 
    228 !!          ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    229 !!          ikbum1 = MAX( ikbu-1, 1 ) 
    230 !!          ikbvm1 = MAX( ikbv-1, 1 ) 
    231 !!          ikbu = MIN( mbathy(ji,jj), mbathy(ji-1,jj) ) 
    232 !!          ikbv = MIN( mbathy(ji,jj), mbathy(ji,jj-1) ) 
    233 !!          ikbumm1 = MAX( ikbu-1, 1 ) 
    234 !!          ikbvmm1 = MAX( ikbv-1, 1 ) 
    235 !!          ikbt = MAX( mbathy(ji,jj), 1 ) 
    236 !!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,ikbumm1) + & 
    237 !!                 bfrua(ji  ,jj) * ub(ji  ,jj  ,ikbum1 ) 
    238 !!          zty2 = bfrva(ji,jj  ) * vb(ji  ,jj  ,ikbvm1) + & 
    239 !!                 bfrva(ji,jj-1) * vb(ji  ,jj-1,ikbvmm1 ) 
     227!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 
     228!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) ) 
     229!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + & 
     230!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 
    240231!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 
    241 !!          en (ji,jj,ikbt) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 
     232!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 
    242233!!       END DO 
    243234!!    END DO 
    … …  
    255246         !                        !* finite Langmuir Circulation depth 
    256247         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 
    257          imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point mbathy 
     248         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land) 
    258249         DO jk = jpkm1, 2, -1 
    259250            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us  
    … …  
    509500               DO ji = fs_2, fs_jpim1   ! vector opt. 
    510501                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   & 
    511                   &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) ) 
     502                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    512503                  zmxlm(ji,jj,jk) = zemxl 
    513504                  zmxld(ji,jj,jk) = zemxl 

  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r2287 r2450  
    371371      DO jj = 1, jpj                ! part independent of the level 
    372372         DO ji = 1, jpi 
    373             zhdep(ji,jj) = fsdepw(ji,jj,mbathy(ji,jj))       ! depth of the ocean 
     373            zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
    374374            zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) ) 
    375375            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj) 
Note: See TracChangeset for help on using the changeset viewer.