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 5972 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2015-12-02T09:52:20+01:00 (8 years ago)
Author:
timgraham
Message:

Upgraded to head of trunk (r5936)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r5967 r5972  
    2828   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    2929   !!---------------------------------------------------------------------- 
    30 #if defined key_zdftke   ||   defined key_esopa 
     30#if defined key_zdftke 
    3131   !!---------------------------------------------------------------------- 
    3232   !!   'key_zdftke'                                   TKE vertical physics 
     
    5353   USE timing         ! Timing 
    5454   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     55#if defined key_agrif 
     56   USE agrif_opa_interp 
     57   USE agrif_opa_update 
     58#endif 
    5559 
    5660   IMPLICIT NONE 
     
    8589   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3) 
    8690 
    87    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2] 
    8891   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    8992   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
    90    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz 
     93   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation 
    9294#if defined key_c1d 
    9395   !                                                                        !!** 1D cfg only  **   ('key_c1d') 
     
    100102#  include "vectopt_loop_substitute.h90" 
    101103   !!---------------------------------------------------------------------- 
    102    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     104   !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
    103105   !! $Id$ 
    104106   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    115117         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          & 
    116118#endif 
    117          &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     &  
    118          &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          & 
    119          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      ) 
     119         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     &  
     120         &      apdlr(jpi,jpj,jpk) ,                                           STAT= zdf_tke_alloc      ) 
    120121         ! 
    121122      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc ) 
     
    173174      !!---------------------------------------------------------------------- 
    174175      ! 
     176#if defined key_agrif  
     177      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) 
     178      IF( .NOT.Agrif_Root() )   CALL Agrif_Tke 
     179#endif 
     180      ! 
    175181      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    176182         avt (:,:,:) = avt_k (:,:,:)  
     
    189195      avmv_k(:,:,:) = avmv(:,:,:)  
    190196      ! 
    191    END SUBROUTINE zdf_tke 
     197#if defined key_agrif 
     198      ! Update child grid f => parent grid  
     199      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only 
     200#endif       
     201     !  
     202  END SUBROUTINE zdf_tke 
    192203 
    193204 
     
    221232      REAL(wp) ::   zzd_up, zzd_lw                  !    -         - 
    222233!!bfr      REAL(wp) ::   zebot                           !    -         - 
    223       INTEGER , POINTER, DIMENSION(:,:  ) :: imlc 
    224       REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc 
    225       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 
     234      INTEGER , POINTER, DIMENSION(:,:  ) ::   imlc 
     235      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc 
     236      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv 
     237      REAL(wp)                            ::   zri  !   local Richardson number 
    226238      !!-------------------------------------------------------------------- 
    227239      ! 
    228240      IF( nn_timing == 1 )  CALL timing_start('tke_tke') 
    229241      ! 
    230       CALL wrk_alloc( jpi,jpj, imlc )    ! integer 
    231       CALL wrk_alloc( jpi,jpj, zhlc )  
    232       CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
     242      CALL wrk_alloc( jpi,jpj,       imlc )    ! integer 
     243      CALL wrk_alloc( jpi,jpj,       zhlc )  
     244      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )  
    233245      ! 
    234246      zbbrau = rn_ebb / rau0       ! Local constant initialisation 
     
    244256         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin 
    245257            DO ji = fs_2, fs_jpim1   ! vector opt. 
    246                en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 
     258               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1) 
    247259            END DO 
    248260         END DO 
     
    265277      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    266278      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
    267 !CDIR NOVERRCHK 
    268279!!    DO jj = 2, jpjm1 
    269 !CDIR NOVERRCHK 
    270280!!       DO ji = fs_2, fs_jpim1   ! vector opt. 
    271281!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 
     
    306316         END DO 
    307317         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    308 !CDIR NOVERRCHK 
    309318         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    310 !CDIR NOVERRCHK 
    311             DO jj = 2, jpjm1 
    312 !CDIR NOVERRCHK 
     319            DO jj = 2, jpjm1 
    313320               DO ji = fs_2, fs_jpim1   ! vector opt. 
    314321                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
     
    332339      ! 
    333340      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form) 
    334          DO jj = 1, jpj                 ! here avmu, avmv used as workspace 
    335             DO ji = 1, jpi 
    336                avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    337                   &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
    338                   &                            / (  fse3uw_n(ji,jj,jk)               & 
    339                   &                              *  fse3uw_b(ji,jj,jk)  ) 
    340                avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   & 
    341                   &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   & 
    342                   &                            / (  fse3vw_n(ji,jj,jk)               & 
    343                   &                              *  fse3vw_b(ji,jj,jk)  ) 
    344             END DO 
    345          END DO 
    346       END DO 
    347       ! 
     341         DO jj = 1, jpjm1 
     342            DO ji = 1, fs_jpim1   ! vector opt. 
     343               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   & 
     344                  &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   & 
     345                  &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) & 
     346                  &                 / (  fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) ) 
     347               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   & 
     348                  &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   & 
     349                  &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) & 
     350                  &                 / (  fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) ) 
     351            END DO 
     352         END DO 
     353      END DO 
     354      ! 
     355      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: compute apdlr 
     356         ! Note that zesh2 is also computed in the next loop. 
     357         ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 
     358         DO jk = 2, jpkm1 
     359            DO jj = 2, jpjm1 
     360               DO ji = fs_2, fs_jpim1   ! vector opt. 
     361                  !                                          ! shear prod. at w-point weightened by mask 
     362                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
     363                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
     364                  !                                          ! local Richardson number 
     365                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 
     366                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  ) 
     367                   
     368               END DO 
     369            END DO 
     370         END DO 
     371         ! 
     372      ENDIF 
     373      !          
    348374      DO jk = 2, jpkm1           !* Matrix and right hand side in en 
    349375         DO jj = 2, jpjm1 
     
    354380               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal 
    355381                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
    356                   !                                                           ! shear prod. at w-point weightened by mask 
    357                zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
    358                   &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
    359                   ! 
     382               !                                   ! shear prod. at w-point weightened by mask 
     383               zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   & 
     384                  &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )     
     385               ! 
    360386               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw) 
    361387               zd_lw(ji,jj,jk) = zzd_lw 
     
    377403         END DO 
    378404      END DO 
    379       ! 
    380       ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    381       DO jj = 2, jpjm1 
     405      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    382406         DO ji = fs_2, fs_jpim1   ! vector opt. 
    383407            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     
    391415         END DO 
    392416      END DO 
    393       ! 
    394       ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    395       DO jj = 2, jpjm1 
     417      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    396418         DO ji = fs_2, fs_jpim1   ! vector opt. 
    397419            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
     
    434456         END DO 
    435457      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability) 
    436 !CDIR NOVERRCHK 
    437458         DO jk = 2, jpkm1 
    438 !CDIR NOVERRCHK 
    439             DO jj = 2, jpjm1 
    440 !CDIR NOVERRCHK 
     459            DO jj = 2, jpjm1 
    441460               DO ji = fs_2, fs_jpim1   ! vector opt. 
    442461                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
     
    453472      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    454473      ! 
    455       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
    456       CALL wrk_dealloc( jpi,jpj, zhlc )  
    457       CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw )  
     474      CALL wrk_dealloc( jpi,jpj,       imlc )    ! integer 
     475      CALL wrk_dealloc( jpi,jpj,       zhlc )  
     476      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )  
    458477      ! 
    459478      IF( nn_timing == 1 )  CALL timing_stop('tke_tke') 
     
    499518      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    500519      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars 
    501       REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      - 
     520      REAL(wp) ::   zdku, zri, zsqen            !   -      - 
    502521      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      - 
    503522      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld 
     
    529548      ENDIF 
    530549      ! 
    531 !CDIR NOVERRCHK 
    532550      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    533 !CDIR NOVERRCHK 
    534          DO jj = 2, jpjm1 
    535 !CDIR NOVERRCHK 
     551         DO jj = 2, jpjm1 
    536552            DO ji = fs_2, fs_jpim1   ! vector opt. 
    537553               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    538                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 
     554               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    539555            END DO 
    540556         END DO 
     
    543559      !                     !* Physical limits for the mixing length 
    544560      ! 
    545       zmxld(:,:,1 ) = zmxlm(:,:,1)   ! surface set to the minimum value  
     561      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value  
    546562      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
    547563      ! 
    548564      SELECT CASE ( nn_mxl ) 
    549565      ! 
     566 !!gm Not sure of that coding for ISF.... 
    550567      ! where wmask = 0 set zmxlm == fse3w 
    551568      CASE ( 0 )           ! bounded by the distance to surface and bottom 
     
    606623            END DO 
    607624         END DO 
    608 !CDIR NOVERRCHK 
    609625         DO jk = 2, jpkm1 
    610 !CDIR NOVERRCHK 
    611             DO jj = 2, jpjm1 
    612 !CDIR NOVERRCHK 
     626            DO jj = 2, jpjm1 
    613627               DO ji = fs_2, fs_jpim1   ! vector opt. 
    614628                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) ) 
     
    630644      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt) 
    631645      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    632 !CDIR NOVERRCHK 
    633646      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points 
    634 !CDIR NOVERRCHK 
    635          DO jj = 2, jpjm1 
    636 !CDIR NOVERRCHK 
     647         DO jj = 2, jpjm1 
    637648            DO ji = fs_2, fs_jpim1   ! vector opt. 
    638649               zsqen = SQRT( en(ji,jj,jk) ) 
     
    660671            DO jj = 2, jpjm1 
    661672               DO ji = fs_2, fs_jpim1   ! vector opt. 
    662                   zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    663                   !                                          ! shear 
    664                   zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   & 
    665                     &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) ) 
    666                   zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   & 
    667                     &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) ) 
    668                   !                                          ! local Richardson number 
    669                   zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear ) 
    670                   zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  ) 
    671 !!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!) 
    672 !!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  ) 
    673                   avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
     673                  avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    674674# if defined key_c1d 
    675                   e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
    676                   e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri 
     675                  e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number 
     676!!gm bug NO zri here.... 
     677!!gm remove the specific diag for c1d ! 
     678                  e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk)                            ! c1d config. : save Ri 
    677679# endif 
    678680              END DO 
Note: See TracChangeset for help on using the changeset viewer.