Changeset 8178


Ignore:
Timestamp:
2017-06-15T08:44:30+02:00 (3 years ago)
Author:
gm
Message:

#1883 (HPC-09) - step-9: final step for the removal of avmu, avmv from the code

Location:
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8143 r8178  
    2525   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields 
    2626   !!---------------------------------------------------------------------- 
    27    USE oce             ! ocean dynamics and tracers  
    28    USE dom_oce         ! ocean space and time domain 
    29    USE dynadv, ONLY: ln_dynadv_vec 
    30    USE zdf_oce         ! ocean vertical physics 
    31    USE zdfdrg          ! ocean vertical physics: top/bottom friction 
    32    USE ldftra          ! lateral physics: eddy diffusivity coef. 
    33    USE ldfdyn          ! lateral physics: eddy viscosity   coef. 
    34    USE sbc_oce         ! Surface boundary condition: ocean fields 
    35    USE sbc_ice         ! Surface boundary condition: ice fields 
    36    USE icb_oce         ! Icebergs 
    37    USE icbdia          ! Iceberg budgets 
    38    USE sbcssr          ! restoring term toward SST/SSS climatology 
    39    USE phycst          ! physical constants 
    40    USE zdfmxl          ! mixed layer 
    41    USE dianam          ! build name of file (routine) 
    42 !   USE zdfddm          ! vertical  physics: double diffusion 
    43    USE diahth          ! thermocline diagnostics 
    44    USE wet_dry         ! wetting and drying 
    45    USE sbcwave         ! wave parameters 
     27   USE oce            ! ocean dynamics and tracers  
     28   USE dom_oce        ! ocean space and time domain 
     29   USE phycst         ! physical constants 
     30   USE dianam         ! build name of file (routine) 
     31   USE diahth         ! thermocline diagnostics 
     32   USE dynadv   , ONLY: ln_dynadv_vec 
     33   USE icb_oce        ! Icebergs 
     34   USE icbdia         ! Iceberg budgets 
     35   USE ldftra         ! lateral physics: eddy diffusivity coef. 
     36   USE ldfdyn         ! lateral physics: eddy viscosity   coef. 
     37   USE sbc_oce        ! Surface boundary condition: ocean fields 
     38   USE sbc_ice        ! Surface boundary condition: ice fields 
     39   USE sbcssr         ! restoring term toward SST/SSS climatology 
     40   USE sbcwave        ! wave parameters 
     41   USE wet_dry        ! wetting and drying 
     42   USE zdf_oce        ! ocean vertical physics 
     43   USE zdfdrg         ! ocean vertical physics: top/bottom friction 
     44   USE zdfmxl         ! mixed layer 
    4645   ! 
    47    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    48    USE in_out_manager  ! I/O manager 
    49    USE diatmb          ! Top,middle,bottom output 
    50    USE dia25h          ! 25h Mean output 
    51    USE iom 
    52    USE ioipsl 
     46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     47   USE in_out_manager ! I/O manager 
     48   USE diatmb         ! Top,middle,bottom output 
     49   USE dia25h         ! 25h Mean output 
     50   USE iom            !  
     51   USE ioipsl         !  
    5352 
    5453#if defined key_lim2 
     
    194193      ENDIF 
    195194          
    196       CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current 
    197       CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current 
     195      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current 
     196      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current 
    198197      IF ( iom_use("sbu") ) THEN 
    199198         DO jj = 1, jpj 
     
    206205      ENDIF 
    207206       
    208       CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current 
    209       CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current 
     207      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current 
     208      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current 
    210209      IF ( iom_use("sbv") ) THEN 
    211210         DO jj = 1, jpj 
     
    229228      ENDIF 
    230229 
    231       CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
    232       CALL iom_put( "avs" , avs                        )    ! S vert. eddy diff. coef. 
    233       CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
     230      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef. 
     231      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef. 
     232      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef. 
    234233 
    235234      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) ) 
     
    247246         END DO 
    248247         CALL lbc_lnk( z2d, 'T', 1. ) 
    249          CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient 
     248         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
    250249         z2d(:,:) = SQRT( z2d(:,:) ) 
    251          CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient 
     250         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient 
    252251      ENDIF 
    253252          
     
    262261            END DO 
    263262         END DO 
    264          CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2) 
     263         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    265264      ENDIF 
    266265 
     
    274273            END DO 
    275274         END DO 
    276          CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2) 
     275         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    277276      ENDIF 
    278277      ! 
    279278      IF ( iom_use("eken") ) THEN 
    280          z3d(:,:,jk) = 0._wp                               !      kinetic energy  
     279         z3d(:,:,jk) = 0._wp  
    281280         DO jk = 1, jpkm1 
    282281            DO jj = 2, jpjm1 
    283282               DO ji = fs_2, fs_jpim1   ! vector opt. 
    284                   zztmp  = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    285                   z3d(ji,jj,jk) = 0.25_wp * zztmp * (                                    & 
    286                      &            un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
    287                      &          + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
    288                      &          + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
    289                      &          + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
     283                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     284                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   & 
     285                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   & 
     286                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   & 
     287                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   ) 
    290288               END DO 
    291289            END DO 
    292290         END DO 
    293291         CALL lbc_lnk( z3d, 'T', 1. ) 
    294          CALL iom_put( "eken", z3d )            
     292         CALL iom_put( "eken", z3d )                 ! kinetic energy 
    295293      ENDIF 
    296294      ! 
     
    304302            z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    305303         END DO 
    306          CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction 
    307          CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum 
     304         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction 
     305         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum 
    308306      ENDIF 
    309307       
    310308      IF( iom_use("u_heattr") ) THEN 
    311          z2d(:,:) = 0.e0  
     309         z2d(:,:) = 0._wp  
    312310         DO jk = 1, jpkm1 
    313311            DO jj = 2, jpjm1 
     
    318316         END DO 
    319317         CALL lbc_lnk( z2d, 'U', -1. ) 
    320          CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction 
     318         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction 
    321319      ENDIF 
    322320 
     
    331329         END DO 
    332330         CALL lbc_lnk( z2d, 'U', -1. ) 
    333          CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction 
     331         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction 
    334332      ENDIF 
    335333 
     
    340338            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 
    341339         END DO 
    342          CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction 
     340         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction 
    343341      ENDIF 
    344342       
     
    353351         END DO 
    354352         CALL lbc_lnk( z2d, 'V', -1. ) 
    355          CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction 
     353         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction 
    356354      ENDIF 
    357355 
    358356      IF( iom_use("v_salttr") ) THEN 
    359          z2d(:,:) = 0.e0  
     357         z2d(:,:) = 0._wp  
    360358         DO jk = 1, jpkm1 
    361359            DO jj = 2, jpjm1 
     
    366364         END DO 
    367365         CALL lbc_lnk( z2d, 'V', -1. ) 
    368          CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction 
    369       ENDIF 
    370  
    371       ! Vertical integral of temperature 
     366         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction 
     367      ENDIF 
     368 
    372369      IF( iom_use("tosmint") ) THEN 
    373          z2d(:,:)=0._wp 
     370         z2d(:,:) = 0._wp 
    374371         DO jk = 1, jpkm1 
    375372            DO jj = 2, jpjm1 
    376373               DO ji = fs_2, fs_jpim1   ! vector opt. 
    377                   z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
     374                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
    378375               END DO 
    379376            END DO 
    380377         END DO 
    381378         CALL lbc_lnk( z2d, 'T', -1. ) 
    382          CALL iom_put( "tosmint", z2d )  
    383       ENDIF 
    384  
    385       ! Vertical integral of salinity 
     379         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature 
     380      ENDIF 
    386381      IF( iom_use("somint") ) THEN 
    387382         z2d(:,:)=0._wp 
     
    389384            DO jj = 2, jpjm1 
    390385               DO ji = fs_2, fs_jpim1   ! vector opt. 
    391                   z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     386                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
    392387               END DO 
    393388            END DO 
    394389         END DO 
    395390         CALL lbc_lnk( z2d, 'T', -1. ) 
    396          CALL iom_put( "somint", z2d )  
    397       ENDIF 
    398  
    399       CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2) 
    400       ! 
    401  
    402       IF (ln_diatmb) THEN      ! If we want tmb values  
    403          CALL dia_tmb  
    404       ENDIF  
    405       IF (ln_dia25h) THEN 
    406          CALL dia_25h( kt ) 
    407       ENDIF  
     391         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity 
     392      ENDIF 
     393 
     394      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2) 
     395      ! 
     396 
     397      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values  
     398           
     399      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging 
    408400 
    409401      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 
     
    729721         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt 
    730722            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    731          CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu 
     723         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm 
    732724            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    733725 
     
    856848      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
    857849      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    858       CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     850      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
    859851      IF( ln_zdfddm ) THEN 
    860852         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef. 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r7753 r8178  
    3737 
    3838   !                      ! Parameter to control the type of lateral viscous operator 
    39    INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in setting the operator 
    40    INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral viscous trend) 
     39   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   !: error in setting the operator 
     40   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   !: without operator (i.e. no lateral viscous trend) 
    4141   !                          !!      laplacian     !    bilaplacian    ! 
    42    INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator 
    43    INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       ! iso-neutral or geopotential operator 
     42   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator 
     43   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator 
    4444 
    45    INTEGER ::   nldf   ! type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
     45   INTEGER, PUBLIC ::   nldf   !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 
    4646 
    4747   !! * Substitutions 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r6140 r8178  
    3737   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90 
    3838 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity 
     40    
    3941   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)  
    4042   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      - 
     
    5355      !!                  ***  ROUTINE dyn_ldf_iso_alloc  *** 
    5456      !!---------------------------------------------------------------------- 
    55       ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     &  
    56          &      zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 
     57      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     &  
     58         &      akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 
    5759         ! 
    5860      IF( dyn_ldf_iso_alloc /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.') 
     
    99101      !! 
    100102      !! ** Action : 
    101       !!        Update (ua,va) arrays with the before geopotential biharmonic 
    102       !!      mixing trend. 
    103       !!        Update (avmu,avmv) to accompt for the diagonal vertical component 
    104       !!      of the rotated operator in dynzdf module 
     103      !!       -(ua,va) updated with the before geopotential harmonic mixing trend 
     104      !!       -(akzu,akzv) to accompt for the diagonal vertical component 
     105      !!                    of the rotated operator in dynzdf module 
    105106      !!---------------------------------------------------------------------- 
    106107      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     
    144145         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
    145146         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
    146   
    147 !!bug 
    148          IF( kt == nit000 ) then 
    149             IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  & 
    150                &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj)) 
    151          endif 
    152 !!end 
    153       ENDIF 
     147         ! 
     148       ENDIF 
    154149 
    155150      !                                                ! =============== 
     
    365360                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     & 
    366361                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) 
    367                ! update avmu (add isopycnal vertical coefficient to avmu) 
    368                ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
    369                avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 
     362               ! vertical mixing coefficient (akzu) 
     363               ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
     364               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 
    370365            END DO 
    371366         END DO 
     
    391386                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
    392387                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
    393                ! update avmv (add isopycnal vertical coefficient to avmv) 
    394                ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
    395                avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 
     388               ! vertical mixing coefficient (akzv) 
     389               ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
     390               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 
    396391            END DO 
    397392         END DO 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r8160 r8178  
    66   !! History :  1.0  !  2005-11  (G. Madec)  Original code 
    77   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
    8    !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option 
    9    !!---------------------------------------------------------------------- 
    10  
    11    !!---------------------------------------------------------------------- 
    12    !!   dyn_zdf       : Update the momentum trend with the vertical diffusion 
     8   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option + avm at t-point 
     9   !!---------------------------------------------------------------------- 
     10 
     11   !!---------------------------------------------------------------------- 
     12   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing 
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce            ! ocean dynamics and tracers variables 
     
    1818   USE zdf_oce        ! ocean vertical physics variables 
    1919   USE zdfdrg         ! vertical physics: top/bottom drag coef. 
    20    USE dynadv   , ONLY: ln_dynadv_vec ! dynamics: advection form 
     20   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form 
     21   USE dynldf    ,ONLY: nldf, np_lap_i   ! dynamics: type of lateral mixing  
     22   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing  
    2123   USE ldfdyn         ! lateral diffusion: eddy viscosity coef. 
    2224   USE trd_oce        ! trends: ocean variables 
     
    2628   USE lib_mpp        ! MPP library 
    2729   USE prtctl         ! Print control 
    28    USE wrk_nemo       ! Memory Allocation 
    2930   USE timing         ! Timing 
    3031 
     
    4950      !!                  ***  ROUTINE dyn_zdf  *** 
    5051      !! 
    51       !! ** Purpose :   compute the vertical ocean dynamics physics. 
    52       !!--------------------------------------------------------------------- 
    53       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    54       ! 
    55       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    56       !!--------------------------------------------------------------------- 
    57       ! 
    58       IF( nn_timing == 1 )   CALL timing_start('dyn_zdf') 
    59       ! 
    60       !                             !* set time step 
    61       IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping) 
    62       ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog) 
    63       ENDIF 
    64  
    65       IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
    66          CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    67          ztrdu(:,:,:) = ua(:,:,:) 
    68          ztrdv(:,:,:) = va(:,:,:) 
    69       ENDIF 
    70       !                             !* compute lateral mixing trend and add it to the general trend 
    71       ! 
    72       CALL dyn_zdf_imp( kt, r2dt ) 
    73       ! 
    74       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    75          ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
    76          ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
    77          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
    78          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    79       ENDIF 
    80       !                                          ! print mean trends (used for debugging) 
    81       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               & 
    82          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    83          ! 
    84       IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf') 
    85       ! 
    86    END SUBROUTINE dyn_zdf 
    87  
    88  
    89    SUBROUTINE dyn_zdf_imp( kt, p2dt ) 
    90       !!---------------------------------------------------------------------- 
    91       !!                  ***  ROUTINE dyn_zdf_imp  *** 
    92       !!                    
    93       !! ** Purpose :   Compute the trend due to the vert. momentum diffusion 
     52      !! ** Purpose :   compute the trend due to the vert. momentum diffusion 
    9453      !!              together with the Leap-Frog time stepping using an  
    9554      !!              implicit scheme. 
     
    10059      !!               - update the after velocity with the implicit vertical mixing. 
    10160      !!      This requires to solver the following system:  
    102       !!         ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ] 
     61      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 
    10362      !!      with the following surface/top/bottom boundary condition: 
    10463      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
    10564      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 
    10665      !! 
    107       !! ** Action :   (ua,va) after velocity  
     66      !! ** Action :   (ua,va)   after velocity  
    10867      !!--------------------------------------------------------------------- 
    10968      INTEGER , INTENT(in) ::  kt     ! ocean time-step index 
    110       REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    111       ! 
    112       INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    113       INTEGER  ::   ikbu, ikbv    ! local integers 
    114       REAL(wp) ::   zzwi, ze3ua   ! local scalars 
    115       REAL(wp) ::   zzws, ze3va   !   -      - 
    116       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwd, zws 
    117       !!---------------------------------------------------------------------- 
    118       ! 
    119       IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp') 
    120       ! 
    121       IF( kt == nit000 ) THEN 
     69      ! 
     70      INTEGER  ::   ji, jj, jk         ! dummy loop indices 
     71      INTEGER  ::   iku, ikv           ! local integers 
     72      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars 
     73      REAL(wp) ::   zzws, ze3va        !   -      - 
     74      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace  
     75      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      - 
     76      !!--------------------------------------------------------------------- 
     77      ! 
     78      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf') 
     79      ! 
     80      IF( kt == nit000 ) THEN       !* initialization 
    12281         IF(lwp) WRITE(numout,*) 
    12382         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
     
    12887         ENDIF 
    12988      ENDIF 
    130       ! 
    131       !              !==  Time step dynamics  ==! 
    132       ! 
    133       IF( ln_dynadv_vec .OR. ln_linssh ) THEN      ! applied on velocity 
    134          DO jk = 1, jpkm1 
    135             ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    136             va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    137          END DO 
    138       ELSE                                         ! applied on thickness weighted velocity 
     89      !                             !* set time step 
     90      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping) 
     91      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog) 
     92      ENDIF 
     93 
     94      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends 
     95         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) )  
     96         ztrdu(:,:,:) = ua(:,:,:) 
     97         ztrdv(:,:,:) = va(:,:,:) 
     98      ENDIF 
     99      ! 
     100      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va) 
     101      ! 
     102      !                    ! time stepping except vertical diffusion 
     103      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity 
     104         DO jk = 1, jpkm1 
     105            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     106            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     107         END DO 
     108      ELSE                                      ! applied on thickness weighted velocity 
    139109         DO jk = 1, jpkm1 
    140110            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
    141                &          + p2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
     111               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
    142112            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
    143                &          + p2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
    144          END DO 
    145       ENDIF 
    146       ! 
    147       !              !==  Apply semi-implicit bottom friction  ==! 
    148       ! 
    149       ! Only needed for semi-implicit bottom friction setup. The explicit 
    150       ! bottom friction has been included in "u(v)a" which act as the R.H.S 
    151       ! column vector of the tri-diagonal matrix equation 
    152       ! 
    153       IF( ln_drgimp ) THEN 
    154          DO jj = 2, jpjm1 
    155             DO ji = 2, jpim1 
    156                ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
    157                ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
    158                avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1) 
    159                avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1) 
    160             END DO 
    161          END DO 
    162          IF ( ln_isfcav ) THEN 
    163             DO jj = 2, jpjm1 
    164                DO ji = 2, jpim1 
    165                   ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
    166                   ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    167                   ! top Cd is masked (=0 outside cavities) no need of test on mik>=2 
    168                   avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu) 
    169                   avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv) 
    170                END DO 
    171             END DO 
    172          END IF 
    173       ENDIF 
    174       ! 
    175       ! With split-explicit free surface, barotropic stress is treated explicitly 
    176       ! Update velocities at the bottom. 
    177       ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
    178       !            not lead to the effective stress seen over the whole barotropic loop.  
    179       ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 
     113               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
     114         END DO 
     115      ENDIF 
     116      !                    ! add top/bottom friction  
     117      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom. 
     118      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
     119      !                not lead to the effective stress seen over the whole barotropic loop.  
     120      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 
    180121      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
    181122         DO jk = 1, jpkm1        ! remove barotropic velocities 
     
    185126         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
    186127            DO ji = fs_2, fs_jpim1   ! vector opt. 
    187                ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    188                ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    189                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
    190                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
    191                ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    192                va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
     128               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     129               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     130               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
     131               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
     132               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     133               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
    193134            END DO 
    194135         END DO 
     
    196137            DO jj = 2, jpjm1         
    197138               DO ji = fs_2, fs_jpim1   ! vector opt. 
    198                   ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
    199                   ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    200                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
    201                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
    202                   ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    203                   va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     139                  iku = miku(ji,jj)         ! top ocean level at u- and v-points  
     140                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     141                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 
     142                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 
     143                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
     144                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
    204145               END DO 
    205146            END DO 
     
    209150      !              !==  Vertical diffusion on u  ==! 
    210151      ! 
    211       ! Matrix and second member construction 
    212       ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    213       ! non zero value at the ocean bottom depending on the bottom friction used. 
    214       ! 
    215       DO jk = 1, jpkm1        ! Matrix 
    216          DO jj = 2, jpjm1  
    217             DO ji = fs_2, fs_jpim1   ! vector opt. 
    218                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
    219                zzwi = - p2dt * avmu(ji,jj,jk  ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) 
    220                zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 
    221                zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk  ) 
    222                zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 
    223                zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    224             END DO 
    225          END DO 
    226       END DO 
    227       DO jj = 2, jpjm1        ! Surface boundary conditions 
     152      !                    !* Matrix construction 
     153      zdt = r2dt * 0.5 
     154      IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu) 
     155         DO jk = 1, jpkm1 
     156            DO jj = 2, jpjm1  
     157               DO ji = fs_2, fs_jpim1   ! vector opt. 
     158                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
     159                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   & 
     160                     &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     161                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   & 
     162                     &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     163                  zwi(ji,jj,jk) = zzwi 
     164                  zws(ji,jj,jk) = zzws 
     165                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     166               END DO 
     167            END DO 
     168         END DO 
     169      ELSE                          ! standard case 
     170         DO jk = 1, jpkm1 
     171            DO jj = 2, jpjm1  
     172               DO ji = fs_2, fs_jpim1   ! vector opt. 
     173                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
     174                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  ) 
     175                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 
     176                  zwi(ji,jj,jk) = zzwi 
     177                  zws(ji,jj,jk) = zzws 
     178                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     179               END DO 
     180            END DO 
     181         END DO 
     182      ENDIF 
     183      ! 
     184      DO jj = 2, jpjm1     !* Surface boundary conditions 
    228185         DO ji = fs_2, fs_jpim1   ! vector opt. 
    229186            zwi(ji,jj,1) = 0._wp 
     
    231188         END DO 
    232189      END DO 
    233  
     190      ! 
     191      !              !==  Apply semi-implicit bottom friction  ==! 
     192      ! 
     193      !     Only needed for semi-implicit bottom friction setup. The explicit 
     194      !     bottom friction has been included in "u(v)a" which act as the R.H.S 
     195      !     column vector of the tri-diagonal matrix equation 
     196      ! 
     197      IF ( ln_drgimp ) THEN      ! implicit bottom friction 
     198         DO jj = 2, jpjm1 
     199            DO ji = 2, jpim1 
     200               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
     201               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
     202               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 
     203            END DO 
     204         END DO 
     205         IF ( ln_isfcav ) THEN   ! top friction (always implicit) 
     206            DO jj = 2, jpjm1 
     207               DO ji = 2, jpim1 
     208                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed 
     209                  iku = miku(ji,jj)       ! ocean top level at u- and v-points  
     210                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point 
     211                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 
     212               END DO 
     213            END DO 
     214         END IF 
     215      ENDIF 
     216      ! 
    234217      ! Matrix inversion starting from the first level 
    235218      !----------------------------------------------------------------------- 
     
    258241         DO ji = fs_2, fs_jpim1   ! vector opt. 
    259242            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
    260             ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     243            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    261244               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    262245         END DO 
     
    285268      !              !==  Vertical diffusion on v  ==! 
    286269      ! 
    287       ! Matrix and second member construction 
    288       ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    289       ! non zero value at the ocean bottom depending on the bottom friction used 
    290       ! 
    291       DO jk = 1, jpkm1        ! Matrix 
    292          DO jj = 2, jpjm1    
    293             DO ji = fs_2, fs_jpim1   ! vector opt. 
    294                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
    295                zzwi = - p2dt * avmv (ji,jj,jk  ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) 
    296                zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 
    297                zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
    298                zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
    299                zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    300             END DO 
    301          END DO 
    302       END DO 
    303       DO jj = 2, jpjm1        ! Surface boundary conditions 
     270      !                       !* Matrix construction 
     271      zdt = r2dt * 0.5 
     272      IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu) 
     273         DO jk = 1, jpkm1 
     274            DO jj = 2, jpjm1    
     275               DO ji = fs_2, fs_jpim1   ! vector opt. 
     276                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
     277                  zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   & 
     278                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     279                  zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   & 
     280                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     281                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     282                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     283                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     284               END DO 
     285            END DO 
     286         END DO 
     287      ELSE                          ! standard case 
     288         DO jk = 1, jpkm1 
     289            DO jj = 2, jpjm1    
     290               DO ji = fs_2, fs_jpim1   ! vector opt. 
     291                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
     292                  zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  ) 
     293                  zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 
     294                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     295                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     296                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
     297               END DO 
     298            END DO 
     299         END DO 
     300      ENDIF 
     301      ! 
     302      DO jj = 2, jpjm1        !* Surface boundary conditions 
    304303         DO ji = fs_2, fs_jpim1   ! vector opt. 
    305304            zwi(ji,jj,1) = 0._wp 
     
    307306         END DO 
    308307      END DO 
     308      !              !==  Apply semi-implicit top/bottom friction  ==! 
     309      ! 
     310      !     Only needed for semi-implicit bottom friction setup. The explicit 
     311      !     bottom friction has been included in "u(v)a" which act as the R.H.S 
     312      !     column vector of the tri-diagonal matrix equation 
     313      ! 
     314      IF( ln_drgimp ) THEN 
     315         DO jj = 2, jpjm1 
     316            DO ji = 2, jpim1 
     317               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
     318               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
     319               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va            
     320            END DO 
     321         END DO 
     322         IF ( ln_isfcav ) THEN 
     323            DO jj = 2, jpjm1 
     324               DO ji = 2, jpim1 
     325                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     326                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point 
     327                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 
     328               END DO 
     329            END DO 
     330         ENDIF 
     331      ENDIF 
    309332 
    310333      ! Matrix inversion 
     
    334357         DO ji = fs_2, fs_jpim1   ! vector opt.           
    335358            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
    336             va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     359            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    337360               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1)  
    338361         END DO 
     
    359382      END DO 
    360383      ! 
    361       IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp') 
    362       ! 
    363    END SUBROUTINE dyn_zdf_imp 
     384      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
     385         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 
     386         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 
     387         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 
     388         DEALLOCATE( ztrdu, ztrdv )  
     389      ENDIF 
     390      !                                          ! print mean trends (used for debugging) 
     391      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               & 
     392         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     393         ! 
     394      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf') 
     395      ! 
     396   END SUBROUTINE dyn_zdf 
    364397 
    365398   !!============================================================================== 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r8160 r8178  
    66   !! history :  1.0  !  2002-06  (G. Madec)  Original code 
    77   !!            3.2  !  2009-07  (G. Madec)  addition of avm 
    8    !!            4.0  !  2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90) 
     8   !!            4.0  !  2017-05  (G. Madec)  avm and drag coef. defined at t-point 
    99   !!---------------------------------------------------------------------- 
    1010   USE par_oce        ! ocean parameters 
     
    4444 
    4545 
    46    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm, avt, avs  !: vertical mixing coefficient (w-point) [m2/s] 
    47 !!gm 
    48    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts  [m2/s] 
    49 !!gm 
    50    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k, avt_k   !: avm, avt computed by turbulent closure alone 
    51    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2] 
    52    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile 
    53    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt 
     46   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm, avt, avs  !: vertical mixing coefficients (w-point) [m2/s] 
     47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k , avt_k  !: Kz computed by turbulent closure alone [m2/s] 
     48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy          [m2/s2] 
     49   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt      [m2/s] 
     50   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile [-] 
    5451 
    5552   !!---------------------------------------------------------------------- 
     
    6562      !!---------------------------------------------------------------------- 
    6663      ! 
    67       ALLOCATE( avm  (jpi,jpj,jpk) , avt  (jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,   & 
    68          &      avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   &  
    69          &      avmb(jpk) , avtb(jpk) ,  avtb_2d(jpi,jpj) ,                    & 
    70          &      avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk)    , STAT = zdf_oce_alloc ) 
     64      ALLOCATE( avm (jpi,jpj,jpk) , avm_k(jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,   & 
     65         &      avt (jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   &  
     66         &      avmb(jpk)         , avtb(jpk)          , avtb_2d(jpi,jpj) , STAT = zdf_oce_alloc ) 
    7167         ! 
    7268      IF( zdf_oce_alloc /= 0 )   CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays') 
  • branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8160 r8178  
    2929   !!            3.7  !  2014-10  (G. Madec)  LDF simplication  
    3030   !!             -   !  2014-12  (G. Madec) remove KPP scheme 
    31    !!             -   !  2015-11  (J. Chanut) free surface simplification 
     31   !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface) 
     32   !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy) 
    3233   !!---------------------------------------------------------------------- 
    3334 
     
    3637   !!---------------------------------------------------------------------- 
    3738   USE step_oce         ! time stepping definition modules 
    38 !!gm to be removed when removing avmu, avmv 
    39    USE zdf_oce        ! ocean vertical physics variables 
    40 !!gm 
    4139   ! 
    4240   USE iom              ! xIOs server 
     
    4846 
    4947   !!---------------------------------------------------------------------- 
    50    !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
     48   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    5149   !! $Id$ 
    5250   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    131129                         CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 
    132130 
    133  
    134 !!gm  ===>>>   TO BE REMOVED   (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri) 
    135       DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
    136          DO jj = 2, jpjm1 
    137             DO ji = 2, jpim1 
    138                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
    139                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    140             END DO 
    141          END DO 
    142       END DO 
    143       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    144 !!gm end 
    145  
    146  
    147131      !  LATERAL  PHYSICS 
    148132      ! 
Note: See TracChangeset for help on using the changeset viewer.