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 5758 for branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Ignore:
Timestamp:
2015-09-24T08:31:40+02:00 (9 years ago)
Author:
gm
Message:

#1593: LDF-ADV, step II.1: phasing the improvements/simplifications of diffusive trend (see wiki)

Location:
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r4990 r5758  
    1515   USE phycst         ! physical constants 
    1616   USE ldfdyn_oce     ! ocean dynamics lateral physics 
    17    USE ldftra_oce     ! ocean tracers  lateral physics 
     17   USE ldftra         ! ocean tracers  lateral physics 
    1818   USE ldfslp         ! lateral mixing: slopes of mixing orientation 
    1919   USE dynldf_bilapg  ! lateral mixing            (dyn_ldf_bilapg routine) 
    … …  
    7373      CASE ( 1 )    ;   CALL dyn_ldf_iso    ( kt )      ! rotated laplacian (except dk[ dk[.] ] part) 
    7474      CASE ( 2 )    ;   CALL dyn_ldf_bilap  ( kt )      ! iso-level bilaplacian 
    75       CASE ( 3 )    ;   CALL dyn_ldf_bilapg ( kt )      ! s-coord. horizontal bilaplacian 
     75!!gm     CASE ( 3 )    ;   CALL dyn_ldf_bilapg ( kt )      ! s-coord. horizontal bilaplacian 
    7676      CASE ( 4 )                                        ! iso-level laplacian + bilaplacian 
    7777         CALL dyn_ldf_lap    ( kt ) 
    … …  
    7979      CASE ( 5 )                                        ! rotated laplacian + bilaplacian (s-coord) 
    8080         CALL dyn_ldf_iso    ( kt ) 
    81          CALL dyn_ldf_bilapg ( kt ) 
     81!!gm         CALL dyn_ldf_bilapg ( kt ) 
    8282      ! 
    8383      CASE ( -1 )                                       ! esopa: test all possibility with control print 
    … …  
    9191                        CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask,   & 
    9292            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    93                         CALL dyn_ldf_bilapg ( kt ) 
     93!!gm                        CALL dyn_ldf_bilapg ( kt ) 
    9494                        CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask,   & 
    9595            &                         tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    … …  
    215215      IF( ierr == 1 )   CALL ctl_stop( 'iso-level in z-coordinate - partial step, not allowed' ) 
    216216      IF( ierr == 2 )   CALL ctl_stop( 'isoneutral bilaplacian operator does not exist' ) 
    217       IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation 
    218          IF( .NOT.lk_ldfslp )   CALL ctl_stop( 'the rotation of the diffusive tensor require key_ldfslp' ) 
    219       ENDIF 
     217      IF( nldf == 1 .OR. nldf == 3 )   l_ldfslp = .TRUE.    ! the rotation needs slope computation 
    220218 
    221219      IF(lwp) THEN 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90

    r4990 r5758  
    11MODULE dynldf_bilapg 
    2    !!====================================================================== 
    3    !!                       ***  MODULE  dynldf_bilapg  *** 
    4    !! Ocean dynamics:  lateral viscosity trend 
    5    !!====================================================================== 
    6    !! History :  OPA  !  1997-07  (G. Madec)  Original code 
    7    !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module 
    8    !!            2.0  !  2004-08  (C. Talandier) New trends organization 
    9    !!---------------------------------------------------------------------- 
    10 #if defined key_ldfslp   ||   defined key_esopa 
    11    !!---------------------------------------------------------------------- 
    12    !!   'key_ldfslp'                              Rotation of mixing tensor 
    13    !!---------------------------------------------------------------------- 
    14    !!   dyn_ldf_bilapg : update the momentum trend with the horizontal part 
    15    !!                    of the horizontal s-coord. bilaplacian diffusion 
    16    !!   ldfguv         :  
    17    !!---------------------------------------------------------------------- 
    18    USE oce             ! ocean dynamics and tracers 
    19    USE dom_oce         ! ocean space and time domain 
    20    USE ldfdyn_oce      ! ocean dynamics lateral physics 
    21    USE zdf_oce         ! ocean vertical physics 
    22    USE ldfslp          ! iso-neutral slopes available 
    23    USE ldftra_oce, ONLY: ln_traldf_iso 
    24    ! 
    25    USE in_out_manager  ! I/O manager 
    26    USE lib_mpp         ! MPP library 
    27    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    28    USE prtctl          ! Print control 
    29    USE wrk_nemo        ! Memory Allocation 
    30    USE timing          ! Timing 
     2   !!============================================================================== 
     3    
     4    
     5   !!    ====>>>   Empty TO BE REMOVED 
     6    
    317 
    32    IMPLICIT NONE 
    33    PRIVATE 
     8   !!============================================================================== 
    349 
    35    PUBLIC   dyn_ldf_bilapg       ! called by step.F90 
    36  
    37    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zfuw, zfvw , zdiu, zdiv   ! 2D workspace (ldfguv) 
    38    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zdju, zdj1u, zdjv, zdj1v  ! 2D workspace (ldfguv) 
    39  
    40    !! * Substitutions 
    41 #  include "domzgr_substitute.h90" 
    42 #  include "ldfdyn_substitute.h90" 
    43    !!---------------------------------------------------------------------- 
    44    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    45    !! $Id$  
    46    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    47    !!---------------------------------------------------------------------- 
    48 CONTAINS 
    49  
    50    INTEGER FUNCTION dyn_ldf_bilapg_alloc() 
    51       !!---------------------------------------------------------------------- 
    52       !!               ***  ROUTINE dyn_ldf_bilapg_alloc  *** 
    53       !!---------------------------------------------------------------------- 
    54       ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) ,     & 
    55          &      zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc ) 
    56          ! 
    57       IF( dyn_ldf_bilapg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays') 
    58    END FUNCTION dyn_ldf_bilapg_alloc 
    59  
    60  
    61    SUBROUTINE dyn_ldf_bilapg( kt ) 
    62       !!---------------------------------------------------------------------- 
    63       !!                   ***  ROUTINE dyn_ldf_bilapg  *** 
    64       !!                       
    65       !! ** Purpose :   Compute the before trend of the horizontal momentum 
    66       !!      diffusion and add it to the general trend of momentum equation. 
    67       !! 
    68       !! ** Method  :   The lateral momentum diffusive trends is provided by a  
    69       !!      a 4th order operator rotated along geopotential surfaces. It is  
    70       !!      computed using before fields (forward in time) and geopotential 
    71       !!      slopes computed in routine inildf. 
    72       !!         -1- compute the geopotential harmonic operator applied to 
    73       !!      (ub,vb) and multiply it by the eddy diffusivity coefficient 
    74       !!      (done by a call to ldfgpu and ldfgpv routines) The result is in 
    75       !!      (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions 
    76       !!      by call to lbc_lnk. 
    77       !!         -2- applied to (zwk1,zwk2) the geopotential harmonic operator 
    78       !!      by a second call to ldfgpu and ldfgpv routines respectively. The 
    79       !!      result is in (zwk3,zwk4) arrays. 
    80       !!         -3- Add this trend to the general trend (ta,sa): 
    81       !!            (ua,va) = (ua,va) + (zwk3,zwk4) 
    82       !! 
    83       !! ** Action  : - Update (ua,va) arrays with the before geopotential 
    84       !!                biharmonic mixing trend. 
    85       !!---------------------------------------------------------------------- 
    86       INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
    87       ! 
    88       INTEGER ::   ji, jj, jk                 ! dummy loop indices 
    89       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwk1, zwk2, zwk3, zwk4 
    90       !!---------------------------------------------------------------------- 
    91       ! 
    92       IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilapg') 
    93       ! 
    94       CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )  
    95       ! 
    96       IF( kt == nit000 ) THEN 
    97          IF(lwp) WRITE(numout,*) 
    98          IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate' 
    99          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 
    100          !                                      ! allocate dyn_ldf_bilapg arrays 
    101          IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') 
    102       ENDIF 
    103  
    104       ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum 
    105       IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 
    106          ! 
    107          DO jk = 1, jpk         ! set the slopes of iso-level 
    108             DO jj = 2, jpjm1 
    109                DO ji = 2, jpim1 
    110                   uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
    111                   vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
    112                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
    113                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
    114                END DO 
    115             END DO 
    116          END DO 
    117          ! Lateral boundary conditions on the slopes 
    118          CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. ) 
    119          CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
    120   
    121 !!bug 
    122          IF( kt == nit000 ) then 
    123             IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  & 
    124                &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj)) 
    125          endif 
    126 !!end 
    127       ENDIF 
    128  
    129       zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0 
    130       zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0 
    131  
    132       ! Laplacian of (ub,vb) multiplied by ahm 
    133       ! --------------------------------------   
    134       CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb) 
    135       !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) ) 
    136       CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions 
    137  
    138       ! Bilaplacian of (ub,vb) 
    139       ! ----------------------  
    140       CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2)  
    141       !                                         ! (output in (zwk3,zwk4) ) 
    142  
    143       ! Update the momentum trends 
    144       ! -------------------------- 
    145       DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends 
    146          DO jk = 1, jpkm1 
    147             DO ji = 2, jpim1 
    148                ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk) 
    149                va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk) 
    150             END DO 
    151          END DO 
    152       END DO 
    153       ! 
    154       CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )  
    155       ! 
    156       IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilapg') 
    157       ! 
    158    END SUBROUTINE dyn_ldf_bilapg 
    159  
    160  
    161    SUBROUTINE ldfguv( pu, pv, plu, plv, kahm ) 
    162       !!---------------------------------------------------------------------- 
    163       !!                  ***  ROUTINE ldfguv  *** 
    164       !!       
    165       !! ** Purpose :   Apply a geopotential harmonic operator to (pu,pv) 
    166       !!      (defined at u- and v-points) and multiply it by the eddy 
    167       !!      viscosity coefficient (if kahm=1). 
    168       !! 
    169       !! ** Method  :   The harmonic operator rotated along geopotential  
    170       !!      surfaces is applied to (pu,pv) using the slopes of geopotential 
    171       !!      surfaces computed in inildf routine. The result is provided in 
    172       !!      (plu,plv) arrays. It is computed in 2 stepv: 
    173       !! 
    174       !!      First step: horizontal part of the operator. It is computed on 
    175       !!      ==========  pu as follows (idem on pv) 
    176       !!      horizontal fluxes : 
    177       !!         zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ] 
    178       !!         zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ] 
    179       !!      take the horizontal divergence of the fluxes (no divided by 
    180       !!      the volume element : 
    181       !!         plu  = di-1[ zftu ] +  dj-1[ zftv ] 
    182       !! 
    183       !!      Second step: vertical part of the operator. It is computed on 
    184       !!      ===========  pu as follows (idem on pv) 
    185       !!      vertical fluxes : 
    186       !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pu ] 
    187       !!              -     e2t     *       wslpi        di[ mi(mk(pu)) ] 
    188       !!              -     e1t     *       wslpj        dj[ mj(mk(pu)) ] 
    189       !!      take the vertical divergence of the fluxes add it to the hori- 
    190       !!      zontal component, divide the result by the volume element and 
    191       !!      if kahm=1, multiply by the eddy diffusivity coefficient: 
    192       !!         plu  = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] } 
    193       !!      else: 
    194       !!         plu  =  1  / (e1t*e2t*e3t) { plu + dk[ zftw ] } 
    195       !! 
    196       !! ** Action : 
    197       !!        plu, plv        : partial harmonic operator applied to 
    198       !!                          pu and pv (all the components except 
    199       !!                          second order vertical derivative term) 
    200       !!---------------------------------------------------------------------- 
    201       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity  
    202       !                                                               ! 2nd call: ahm x these fields 
    203       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial harmonic operator applied to 
    204       !                                                               ! pu and pv (all the components except 
    205       !                                                               ! second order vertical derivative term) 
    206       INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call 
    207       ! 
    208       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    209       REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar 
    210       REAL(wp) ::   zcoef0, zcoef3, zcoef4               !   -      - 
    211       REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      - 
    212       REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
    213       ! 
    214       REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v 
    215       !!---------------------------------------------------------------------- 
    216       ! 
    217       IF( nn_timing == 1 )  CALL timing_start('ldfguv') 
    218       ! 
    219       CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )  
    220       ! 
    221       !                               ! ********** !   ! =============== 
    222       DO jk = 1, jpkm1                ! First step !   ! Horizontal slab 
    223          !                            ! ********** !   ! =============== 
    224  
    225          ! I.1 Vertical gradient of pu and pv at level jk and jk+1 
    226          ! ------------------------------------------------------- 
    227          ! surface boundary condition: zdku(jk=1)=zdku(jk=2) 
    228          !                             zdkv(jk=1)=zdkv(jk=2) 
    229  
    230          zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1) 
    231          zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1) 
    232  
    233          IF( jk == 1 ) THEN 
    234             zdku(:,:) = zdk1u(:,:) 
    235             zdkv(:,:) = zdk1v(:,:) 
    236          ELSE 
    237             zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk) 
    238             zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk) 
    239          ENDIF 
    240  
    241          !                                -----f----- 
    242          ! I.2 Horizontal fluxes on U          | 
    243          ! ------------------------===     t   u   t 
    244          !                                     | 
    245          ! i-flux at t-point              -----f----- 
    246          DO jj = 1, jpjm1 
    247             DO ji = 2, jpi 
    248                zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
    249  
    250                zmkt  = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    251                               + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    252  
    253                zcof1 = -e2t(ji,jj) * zmkt   & 
    254                      * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    255  
    256                ziut(ji,jj) = tmask(ji,jj,jk) *   & 
    257                            (  zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   & 
    258                             + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
    259                                        +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) 
    260             END DO 
    261          END DO 
    262  
    263          ! j-flux at f-point 
    264          DO jj = 1, jpjm1 
    265             DO ji = 1, jpim1 
    266                zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
    267  
    268                zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
    269                               + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
    270  
    271                zcof2 = -e1f(ji,jj) * zmkf   & 
    272                      * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    273  
    274                zjuf(ji,jj) = fmask(ji,jj,jk) *   & 
    275                            (  zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   & 
    276                             + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     & 
    277                                        +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) 
    278             END DO 
    279          END DO 
    280  
    281          !                                 |   t   | 
    282          ! I.3 Horizontal fluxes on V      |       | 
    283          ! ------------------------===     f---v---f 
    284          !                                 |       | 
    285          ! i-flux at f-point               |   t   | 
    286          DO jj = 1, jpjm1 
    287             DO ji = 1, jpim1 
    288                zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
    289  
    290                zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
    291                               + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    292  
    293                zcof1 = -e2f(ji,jj) * zmkf   & 
    294                      * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    295  
    296                zivf(ji,jj) = fmask(ji,jj,jk) *   & 
    297                            (  zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   & 
    298                             + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     & 
    299                                        +zdk1u(ji,jj) + zdku (ji+1,jj) )  ) 
    300             END DO 
    301          END DO 
    302  
    303          ! j-flux at t-point 
    304          DO jj = 2, jpj 
    305             DO ji = 1, jpim1 
    306                zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
    307  
    308                zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    309                               + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    310  
    311                zcof2 = -e1t(ji,jj) * zmkt   & 
    312                      * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    313  
    314                zjvt(ji,jj) = tmask(ji,jj,jk) *   & 
    315                            (  zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   & 
    316                             + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     & 
    317                                        +zdk1u(ji,jj-1) + zdku (ji,jj) )  ) 
    318             END DO 
    319          END DO 
    320  
    321  
    322          ! I.4 Second derivative (divergence) (not divided by the volume) 
    323          ! --------------------- 
    324  
    325          DO jj = 2, jpjm1 
    326             DO ji = 2, jpim1 
    327                plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj  )   & 
    328                              + zjuf (ji  ,jj) - zjuf (ji,jj-1) 
    329                plv(ji,jj,jk) = zivf (ji,jj  ) - zivf (ji-1,jj)   & 
    330                              + zjvt (ji,jj+1) - zjvt (ji,jj  )  
    331             END DO 
    332          END DO 
    333  
    334          !                                             ! =============== 
    335       END DO                                           !   End of slab 
    336       !                                                ! =============== 
    337  
    338       !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,, 
    339  
    340       !                             ! ************ !   ! =============== 
    341       DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab 
    342          !                          ! ************ !   ! =============== 
    343  
    344          ! II.1 horizontal (pu,pv) gradients 
    345          ! --------------------------------- 
    346  
    347          DO jk = 1, jpk 
    348             DO ji = 2, jpi 
    349                ! i-gradient of u at jj 
    350                zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) ) 
    351                ! j-gradient of u and v at jj 
    352                zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) ) 
    353                zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) ) 
    354                ! j-gradient of u and v at jj+1 
    355                zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) ) 
    356                zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) ) 
    357             END DO 
    358          END DO 
    359          DO jk = 1, jpk 
    360             DO ji = 1, jpim1 
    361                ! i-gradient of v at jj 
    362                zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) ) 
    363             END DO 
    364          END DO 
    365  
    366  
    367          ! II.2 Vertical fluxes 
    368          ! -------------------- 
    369  
    370          ! Surface and bottom vertical fluxes set to zero 
    371  
    372          zfuw(:, 1 ) = 0.e0 
    373          zfvw(:, 1 ) = 0.e0 
    374          zfuw(:,jpk) = 0.e0 
    375          zfvw(:,jpk) = 0.e0 
    376  
    377          ! interior (2=<jk=<jpk-1) on pu field 
    378  
    379          DO jk = 2, jpkm1 
    380             DO ji = 2, jpim1 
    381                ! i- and j-slopes at uw-point 
    382                zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
    383                zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 
    384                ! coef. for the vertical dirative 
    385                zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk)   & 
    386                       * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) 
    387                ! weights for the i-k, j-k averaging at t- and f-points, resp. 
    388                zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   & 
    389                              + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. ) 
    390                zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   & 
    391                              + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. ) 
    392                ! coef. for the horitontal derivative 
    393                zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 
    394                zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj 
    395                ! vertical flux on u field 
    396                zfuw(ji,jk) = umask(ji,jj,jk) *   & 
    397                            (  zcoef0 * ( pu  (ji,jj,jk-1) - pu  (ji,jj,jk) )   & 
    398                             + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     & 
    399                                         +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   & 
    400                             + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     & 
    401                                         +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) ) 
    402             END DO 
    403          END DO 
    404  
    405          ! interior (2=<jk=<jpk-1) on pv field 
    406  
    407          DO jk = 2, jpkm1 
    408             DO ji = 2, jpim1 
    409                ! i- and j-slopes at vw-point 
    410                zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 
    411                zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) ) 
    412                ! coef. for the vertical derivative 
    413                zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk)   & 
    414                       * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) 
    415                ! weights for the i-k, j-k averaging at f- and t-points, resp. 
    416                zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   & 
    417                              + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. ) 
    418                zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   & 
    419                              + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. ) 
    420                ! coef. for the horizontal derivatives 
    421                zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi 
    422                zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj 
    423                ! vertical flux on pv field 
    424                zfvw(ji,jk) = vmask(ji,jj,jk) *   & 
    425                            (  zcoef0 * ( pv  (ji,jj,jk-1) - pv  (ji,jj,jk) )   & 
    426                             + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     & 
    427                                         +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
    428                             + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
    429                                         +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )  ) 
    430             END DO 
    431          END DO 
    432  
    433  
    434          ! II.3 Divergence of vertical fluxes added to the horizontal divergence 
    435          ! --------------------------------------------------------------------- 
    436          IF( (kahm -nkahm_smag) ==1 ) THEN 
    437             ! multiply the laplacian by the eddy viscosity coefficient 
    438             DO jk = 1, jpkm1 
    439                DO ji = 2, jpim1 
    440                   ! eddy coef. divided by the volume element 
    441                   zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) 
    442                   zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) 
    443                   ! vertical divergence 
    444                   zuav = zfuw(ji,jk) - zfuw(ji,jk+1) 
    445                   zvav = zfvw(ji,jk) - zfvw(ji,jk+1) 
    446                   ! harmonic operator applied to (pu,pv) and multiply by ahm 
    447                   plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur 
    448                   plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr 
    449                END DO 
    450             END DO 
    451          ELSEIF( (kahm +nkahm_smag ) == 2 ) THEN 
    452             ! second call, no multiplication 
    453             DO jk = 1, jpkm1 
    454                DO ji = 2, jpim1 
    455                   ! inverse of the volume element 
    456                   zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) ) 
    457                   zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) ) 
    458                   ! vertical divergence 
    459                   zuav = zfuw(ji,jk) - zfuw(ji,jk+1) 
    460                   zvav = zfvw(ji,jk) - zfvw(ji,jk+1) 
    461                   ! harmonic operator applied to (pu,pv)  
    462                   plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur 
    463                   plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr 
    464                END DO 
    465             END DO 
    466          ELSE 
    467             IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm 
    468             IF(lwp)WRITE(numout,*) '         We stop' 
    469             STOP 'ldfguv' 
    470          ENDIF 
    471          !                                             ! =============== 
    472       END DO                                           !   End of slab 
    473       !                                                ! =============== 
    474  
    475       CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )  
    476       ! 
    477       IF( nn_timing == 1 )  CALL timing_stop('ldfguv') 
    478       ! 
    479    END SUBROUTINE ldfguv 
    480  
    481 #else 
    482    !!---------------------------------------------------------------------- 
    483    !!   Dummy module :                         NO rotation of mixing tensor 
    484    !!---------------------------------------------------------------------- 
    485 CONTAINS 
    486    SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine 
    487       INTEGER, INTENT(in) :: kt 
    488       WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt 
    489    END SUBROUTINE dyn_ldf_bilapg 
    490 #endif 
    491  
    492    !!====================================================================== 
    49310END MODULE dynldf_bilapg 

  • branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r4990 r5758  
    99   !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion 
    1010   !!---------------------------------------------------------------------- 
    11 #if defined key_ldfslp   ||   defined key_esopa 
    12    !!---------------------------------------------------------------------- 
    13    !!   'key_ldfslp'                      slopes of the direction of mixing 
     11 
    1412   !!---------------------------------------------------------------------- 
    1513   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part 
    … …  
    2018   USE dom_oce         ! ocean space and time domain 
    2119   USE ldfdyn_oce      ! ocean dynamics lateral physics 
    22    USE ldftra_oce      ! ocean tracer   lateral physics 
     20   USE ldftra          ! lateral physics: eddy diffusivity & EIV coefficients 
    2321   USE zdf_oce         ! ocean vertical physics 
    2422   USE ldfslp          ! iso-neutral slopes  
    … …  
    106104      !!      of the rotated operator in dynzdf module 
    107105      !!---------------------------------------------------------------------- 
    108       ! 
    109106      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    110107      ! 
    … …  
    189186                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    190187 
    191                   zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     188                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    192189    
    193190                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    … …  
    204201                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    205202 
    206                   zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     203                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    207204 
    208205                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    … …  
    221218                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
    222219 
    223                zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     220               zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    224221 
    225222               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
    … …  
    242239                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    243240 
    244                zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     241               zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    245242 
    246243               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   & 
    … …  
    259256                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    260257 
    261                   zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     258                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    262259 
    263260                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
    … …  
    274271                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    275272 
    276                   zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     273                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    277274 
    278275                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
    … …  
    358355         DO jk = 2, jpkm1 
    359356            DO ji = 2, jpim1 
    360                zcoef0= 0.5 * aht0 * umask(ji,jj,jk) 
    361  
     357               zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk) 
     358               ! 
    362359               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
    363360               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 
    364  
     361               ! 
    365362               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   & 
    366363                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. ) 
    … …  
    376373                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) 
    377374               ! update avmu (add isopycnal vertical coefficient to avmu) 
    378                ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
    379                avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0 
     375               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
     376               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 
    380377            END DO 
    381378         END DO 
    … …  
    384381         DO jk = 2, jpkm1 
    385382            DO ji = 2, jpim1 
    386                zcoef0= 0.5 * aht0 * vmask(ji,jj,jk) 
     383               zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk) 
    387384 
    388385               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 
    … …  
    398395               ! vertical flux on v field 
    399396               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     & 
    400                                        +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
    401                            + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
    402                                        +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
     397                  &                    +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
     398                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
     399                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
    403400               ! update avmv (add isopycnal vertical coefficient to avmv) 
    404                ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
    405                avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0 
     401               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
     402               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 
    406403            END DO 
    407404         END DO 
    … …  
    413410            DO ji = 2, jpim1 
    414411               ! volume elements 
    415                zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    416                zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
     412               zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 
     413               zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 
    417414               ! part of the k-component of isopycnal momentum diffusive trends 
    418415               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 
    … …  
    432429   END SUBROUTINE dyn_ldf_iso 
    433430 
    434 # else 
    435    !!---------------------------------------------------------------------- 
    436    !!   Dummy module                           NO rotation of mixing tensor 
    437    !!---------------------------------------------------------------------- 
    438 CONTAINS 
    439    SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine 
    440       INTEGER, INTENT(in) :: kt 
    441       WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt 
    442    END SUBROUTINE dyn_ldf_iso 
    443 #endif 
    444  
    445431   !!====================================================================== 
    446432END MODULE dynldf_iso 
Note: See TracChangeset for help on using the changeset viewer.