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

Ignore:
Timestamp:
2017-10-16T16:19:11+02:00 (7 years ago)
Author:
gm
Message:

#1963 : Correct implementation of CMIP6 KE dissipation + move EKE<=>PE diag in traadv_eiv (in both cases, no more need of ln_trd_KE=T)

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

Legend:

Unmodified
Added
Removed

  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90

    r4990 r8627  
    1818   USE oce             ! ocean dynamics and tracers 
    1919   USE dom_oce         ! ocean space and time domain 
     20   USE phycst          ! physical constants 
    2021   USE ldfdyn_oce      ! ocean dynamics: lateral physics 
    2122   ! 
    2223   USE in_out_manager  ! I/O manager 
     24   USE iom             ! I/O library 
    2325   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2426   USE wrk_nemo        ! Memory Allocation 
    … …  
    7577      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    7678      ! 
    77       INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    78       REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar 
    79       REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv 
    80       REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv 
     79      INTEGER  ::   ji, jj, jk                       ! dummy loop indices 
     80      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v, zzz   ! local scalar 
     81      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcu, zcv 
     82      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zuf, zut, zlu, zlv 
     83      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z2d   ! 2D workspace  
    8184      !!---------------------------------------------------------------------- 
    8285      ! 
    … …  
    112115            DO jj = 2, jpjm1 
    113116               DO ji = fs_2, fs_jpim1   ! vector opt. 
    114                   zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   & 
    115                      &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj) 
     117                  zlu(ji,jj,jk) = - (   zuf(ji  ,jj,jk) -  zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   & 
     118                     &            + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) /  e1u(ji,jj) 
    116119    
    117                   zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
    118                      &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj) 
     120                  zlv(ji,jj,jk) = + (   zuf(ji,jj  ,jk) -  zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   & 
     121                     &            + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) /  e2v(ji,jj) 
    119122               END DO 
    120123            END DO 
    … …  
    123126               DO ji = fs_2, fs_jpim1   ! vector opt. 
    124127                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   & 
    125                      &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj) 
     128                     &            + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj) 
    126129    
    127130                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   & 
    128                      &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj) 
     131                     &            + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj) 
    129132               END DO   
    130133            END DO   
    … …  
    133136      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions 
    134137 
    135           
     138      IF( iom_use('dispkexyfo') ) THEN   ! ocean kinetic energy dissipation per unit area 
     139         !                               ! due to xy friction (xy=lateral)  
     140         !   see NEMO_book appendix C, §C.7.2 (N.B. here averaged at t-points) 
     141         !   local dissipation of KE at t-point due to bilaplacian operator is given by : 
     142         !      + ahmu mi( zlu**2 ) + mj( ahmv zlv**2 ) 
     143         !   Note that a sign + is used as in v3.6 ahm is negative for bilaplacian viscosity 
     144         ! 
     145         ! NB: ahm is negative when bilaplacian is used 
     146         ALLOCATE( z2d(jpi,jpj) ) 
     147         z2d(:,:) = 0._wp 
     148         DO jk = 1, jpkm1 
     149            DO jj = 2, jpjm1 
     150               DO ji = 2, jpim1 
     151                  z2d(ji,jj) = z2d(ji,jj)                                                                  & 
     152                     &  +  (  fsahmu(ji,jj,jk)*zlu(ji,jj,jk)**2 + fsahmu(ji-1,jj,jk)*zlu(ji-1,jj,jk)**2    & 
     153                     &      + fsahmv(ji,jj,jk)*zlv(ji,jj,jk)**2 + fsahmv(ji,jj-1,jk)*zlv(ji,jj-1,jk)**2  ) * tmask(ji,jj,jk) 
     154               END DO 
     155            END DO 
     156         END DO 
     157         zzz = 0.5_wp * rau0 
     158         z2d(:,:) = zzz * z2d(:,:)  
     159         CALL lbc_lnk( z2d,'T', 1. ) 
     160         CALL iom_put( 'dispkexyfo', z2d ) 
     161         DEALLOCATE( z2d ) 
     162      ENDIF 
     163       
     164    
     165      ! Third derivative 
     166      ! ---------------- 
     167      ! 
    136168      DO jk = 1, jpkm1 
    137     
    138          ! Third derivative 
    139          ! ---------------- 
    140           
     169         ! 
    141170         ! Multiply by the eddy viscosity coef. (at u- and v-points) 
    142          zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 
    143  
    144          zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 
    145           
     171         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)  
     172         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk) 
     173         !          
    146174         ! Contravariant "laplacian" 
    147175         zcu(:,:) = e1u(:,:) * zlu(:,:,jk) 
    … …  
    152180            DO ji = 1, fs_jpim1   ! vector opt. 
    153181               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      & 
    154                   &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   & 
    155                   &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) ) 
     182                  &                               - zcu(ji  ,jj+1) + zcu(ji,jj)  )   & 
     183                  &       * fse3f(ji,jj,jk) * r1_e12f(ji,jj) 
    156184            END DO   
    157185         END DO   
    … …  
    175203      END DO 
    176204 
    177  
    178205      ! boundary conditions on the laplacian curl and div (zuf,zut) 
    179206!!bug gm no need to do this 2 following lbc... 
    … …  
    181208      CALL lbc_lnk( zut, 'T', 1. ) 
    182209 
    183       DO jk = 1, jpkm1       
    184     
    185          ! Bilaplacian 
    186          ! ----------- 
    187  
     210      DO jk = 1, jpkm1               ! Bilaplacian 
    188211         DO jj = 2, jpjm1 
    189212            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    197220                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj) 
    198221               ! add it to the general momentum trends 
    199                ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 
    200                va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 
    201             END DO 
    202          END DO 
    203  
     222               ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
     223               va(ji,jj,jk) = va(ji,jj,jk) + zva 
     224            END DO 
     225         END DO 
    204226         !                                             ! =============== 
    205227      END DO                                           !   End of slab 

  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90

    r4990 r8627  
    1717   USE oce             ! ocean dynamics and tracers 
    1818   USE dom_oce         ! ocean space and time domain 
     19   USE phycst          ! physical constants 
    1920   USE ldfdyn_oce      ! ocean dynamics: lateral physics 
    2021   USE zdf_oce         ! ocean vertical physics 
    2122   ! 
    2223   USE in_out_manager  ! I/O manager 
     24   USE iom             ! I/O library 
    2325   USE timing          ! Timing 
    2426 
    … …  
    6264      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
    6365      REAL(wp) ::   zua, zva, ze2u, ze1v   ! local scalars 
     66      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z2d   ! 2D workspace  
    6467      !!---------------------------------------------------------------------- 
    6568      ! 
    … …  
    9396      END DO                                           !   End of slab 
    9497      !                                                ! =============== 
     98       
     99      IF( iom_use('dispkexyfo') ) THEN   ! ocean Kinetic Energy dissipation per unit area 
     100         !                               ! due to lateral friction (xy=lateral)  
     101         !   see NEMO_book appendix C, §C.7.2 (N.B. here averaged at t-points) 
     102         !   local dissipation of KE at t-point due to laplacian operator is given by : 
     103         !      - ahmt hdivb**2 - mi( mj(ahmf rotb**2 e1f*e2f*e3t) ) / (e1e2t*e3t) 
     104         ! 
     105         ALLOCATE( z2d(jpi,jpj) ) 
     106         z2d(:,:) = 0._wp 
     107         DO jk = 1, jpkm1 
     108            DO jj = 2, jpjm1 
     109               DO ji = 2, jpim1 
     110                  z2d(ji,jj) = z2d(ji,jj)          -   (                                                         & 
     111                     &   hdivb(ji,jj,jk)**2 * fsahmt(ji,jj,jk) * fse3t_n(ji,jj,jk)                               & 
     112                     & + 0.25_wp * (                                                                             & 
     113                     &   rotb (ji  ,jj  ,jk)**2 * fsahmf(ji  ,jj  ,jk) * e12f(ji  ,jj  ) * fse3f(ji  ,jj  ,jk)   & 
     114                     & + rotb (ji-1,jj  ,jk)**2 * fsahmf(ji-1,jj  ,jk) * e12f(ji-1,jj  ) * fse3f(ji-1,jj  ,jk)   & 
     115                     & + rotb (ji  ,jj-1,jk)**2 * fsahmf(ji  ,jj-1,jk) * e12f(ji  ,jj-1) * fse3f(ji  ,jj-1,jk)   & 
     116                     & + rotb (ji-1,jj-1,jk)**2 * fsahmf(ji-1,jj-1,jk) * e12f(ji-1,jj-1) * fse3f(ji-1,jj-1,jk)   & 
     117                     &             ) * r1_e12t(ji,jj)  ) * tmask(ji,jj,jk) 
     118               END DO 
     119            END DO 
     120         END DO 
     121         z2d(:,:) = rau0 * z2d(:,:)  
     122         CALL lbc_lnk( z2d,'T', 1. ) 
     123         CALL iom_put( 'dispkexyfo', z2d ) 
     124         DEALLOCATE( z2d ) 
     125      ENDIF 
     126 
    95127      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap') 
    96128      ! 

  • branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r6751 r8627  
    2020   USE zdf_oce         ! ocean vertical physics 
    2121   USE phycst          ! physical constants 
     22   USE dynadv          ! dynamics: vector invariant versus flux form 
     23   USE dynspg_oce, ONLY: lk_dynspg_ts 
     24   USE zdfbfr          ! Bottom friction setup 
     25   ! 
    2226   USE in_out_manager  ! I/O manager 
     27   USE iom             ! I/O library 
    2328   USE lib_mpp         ! MPP library 
    24    USE zdfbfr          ! Bottom friction setup 
    2529   USE wrk_nemo        ! Memory Allocation 
    2630   USE timing          ! Timing 
    27    USE dynadv          ! dynamics: vector invariant versus flux form 
    28    USE dynspg_oce, ONLY: lk_dynspg_ts 
    2931 
    3032   IMPLICIT NONE 
    … …  
    6971      INTEGER  ::   ikbu, ikbv   ! local integers 
    7072      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
    71       REAL(wp) ::   ze3ua, ze3va 
     73      REAL(wp) ::   ze3ua, ze3va, zzz 
     74      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d 
    7275      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
    7376      !!---------------------------------------------------------------------- 
    … …  
    257260      END DO 
    258261 
    259 #if ! defined key_dynspg_ts 
    260       ! Normalization to obtain the general momentum trend ua 
    261       DO jk = 1, jpkm1 
    262          DO jj = 2, jpjm1    
    263             DO ji = fs_2, fs_jpim1   ! vector opt. 
    264                ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 
    265             END DO 
    266          END DO 
    267       END DO 
    268 #endif 
    269  
    270262      ! 3. Vertical diffusion on v 
    271263      ! --------------------------- 
    … …  
    357349      END DO 
    358350 
    359       ! Normalization to obtain the general momentum trend va 
     351      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area 
     352         !                                ! due to v friction (v=vertical)  
     353         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points) 
     354         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of  
     355         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured). 
     356         CALL wrk_alloc(jpi,jpj,   z2d ) 
     357         z2d(:,:) = 0._wp 
     358         DO jk = 1, jpkm1 
     359            DO jj = 2, jpjm1 
     360               DO ji = 2, jpim1 
     361                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  & 
     362                     &   avmu(ji  ,jj,jk) * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) )**2 / fse3uw(ji  ,jj,jk) * wumask(ji  ,jj,jk)   & 
     363                     & + avmu(ji-1,jj,jk) * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / fse3uw(ji-1,jj,jk) * wumask(ji-1,jj,jk)   & 
     364                     & + avmv(ji,jj  ,jk) * ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) )**2 / fse3vw(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   & 
     365                     & + avmv(ji,jj-1,jk) * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / fse3vw(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   & 
     366                     &                        ) 
     367               END DO 
     368            END DO 
     369         END DO 
     370         zzz= - 0.5_wp* rau0           ! caution sign minus here 
     371         z2d(:,:) = zzz * z2d(:,:)  
     372         CALL lbc_lnk( z2d,'T', 1. ) 
     373         CALL iom_put( 'dispkevfo', z2d ) 
     374         CALL wrk_dealloc(jpi,jpj,   z2d ) 
     375      ENDIF 
     376 
    360377#if ! defined key_dynspg_ts 
     378!!gm this can be removed if tranxt is changed like in the trunk so that implicit outcome with  
     379!!gm the after velocity, not a trend 
     380      ! Normalization to obtain the general momentum trend ua 
    361381      DO jk = 1, jpkm1 
    362382         DO jj = 2, jpjm1    
    363383            DO ji = fs_2, fs_jpim1   ! vector opt. 
     384               ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 
    364385               va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt 
    365386            END DO 
    … …  
    393414      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)  
    394415      ! 
     416      ! 
    395417      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp') 
    396418      ! 
Note: See TracChangeset for help on using the changeset viewer.