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 6383 – NEMO

Changeset 6383


Ignore:
Timestamp:
2016-03-10T17:06:52+01:00 (8 years ago)
Author:
acc
Message:

#1689. 2016WP-SIMPLIF-5. Code changes to ldfdyn.F90 and namelist additions to re-implement Smagorinsky viscosity

Location:
branches/2016/dev_r6381_SIMPLIF_5/NEMOGCM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6381_SIMPLIF_5/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6152 r6383  
    815815   rn_aht_0        = 2000.     !  lateral eddy diffusivity   (lap. operator) [m2/s] 
    816816   rn_bht_0        = 1.e+12    !  lateral eddy diffusivity (bilap. operator) [m4/s] 
     817   ! 
     818                               !  Smagorinsky settings: 
     819   rn_csmc       = 3.5         !  Smagorinsky constant of proportionality 
     820   rn_cfacmin    = 1.0         !  multiplier of theorectical lower limit 
     821   rn_cfacmax    = 1.0         !  multiplier of theorectical upper limit 
    817822/ 
    818823!----------------------------------------------------------------------- 
  • branches/2016/dev_r6381_SIMPLIF_5/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90

    r6140 r6383  
    4242   REAL(wp), PUBLIC ::   rn_ahm_b        !: lateral laplacian background eddy viscosity [m2/s] 
    4343   REAL(wp), PUBLIC ::   rn_bhm_0        !: lateral bilaplacian eddy viscosity          [m4/s] 
     44                                         !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 
     45                                         !! will be computed. 
     46   REAL(wp), PUBLIC ::   rn_csmc         !: Smagorinsky constant of proportionality  
     47   REAL(wp), PUBLIC ::   rn_cfacmin      !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 
     48   REAL(wp), PUBLIC ::   rn_cfacmax      !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 
    4449 
    4550   LOGICAL , PUBLIC ::   l_ldfdyn_time   !: flag for time variation of the lateral eddy viscosity coef. 
    4651 
    4752   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy diffusivity coef. at U- and V-points   [m2/s or m4/s] 
     53   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dtensq       !: horizontal tension squared         (Smagorinsky only) 
     54   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only) 
     55   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2            
    4856 
    4957   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12 
    5058   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4 
     59   REAL(wp) ::   r1_8    = 0.125_wp          ! =1/8 
    5160   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 ) 
    5261 
     
    7988      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator 
    8089      !!                                                           or |u|e^3/12 bilaplacian operator ) 
    81       !!---------------------------------------------------------------------- 
    82       INTEGER  ::   jk                ! dummy loop indices 
     90      !!                  = 32    = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)   
     91      !!                                                           (   L^2|D|      laplacian operator 
     92      !!                                                           or  L^4|D|/8  bilaplacian operator ) 
     93      !!---------------------------------------------------------------------- 
     94      INTEGER  ::   ji, jj, jk        ! dummy loop indices 
    8395      INTEGER  ::   ierr, inum, ios   ! local integer 
    8496      REAL(wp) ::   zah0              ! local scalar 
     
    8698      NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp,                  & 
    8799         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   & 
    88          &                 nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 
     100         &                 nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0,   & 
     101         &                 rn_csmc      , rn_cfacmin, rn_cfacmax 
    89102      !!---------------------------------------------------------------------- 
    90103      ! 
     
    115128         WRITE(numout,*) '      coefficients :' 
    116129         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t 
    117          WRITE(numout,*) '         lateral laplacian eddy viscosity     rn_ahm_0_lap  = ', rn_ahm_0, ' m2/s' 
     130         WRITE(numout,*) '         lateral laplacian eddy viscosity     rn_ahm_0      = ', rn_ahm_0, ' m2/s' 
    118131         WRITE(numout,*) '         background viscosity (iso case)      rn_ahm_b      = ', rn_ahm_b, ' m2/s' 
    119          WRITE(numout,*) '         lateral bilaplacian eddy viscosity   rn_ahm_0_blp  = ', rn_bhm_0, ' m4/s' 
     132         WRITE(numout,*) '         lateral bilaplacian eddy viscosity   rn_bhm_0      = ', rn_bhm_0, ' m4/s' 
     133         WRITE(numout,*) '      smagorinsky settings (nn_ahm_ijk_t  = 32) :' 
     134         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc 
     135         WRITE(numout,*) '         factor multiplier for theorectical lower limit for ' 
     136         WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_cfacmin    = ', rn_cfacmin 
     137         WRITE(numout,*) '         factor multiplier for theorectical lower upper for ' 
     138         WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_cfacmax    = ', rn_cfacmax 
    120139      ENDIF 
    121140 
     
    208227            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90 
    209228            ! 
     229         CASE(  32  )       !==  time varying 3D field  ==! 
     230            IF(lwp) WRITE(numout,*) '          momentum mixing coef. = F( latitude, longitude, depth , time )' 
     231            IF(lwp) WRITE(numout,*) '             proportional to the local deformation rate and gridscale (Smagorinsky)' 
     232            IF(lwp) WRITE(numout,*) '                                                             : L^2|D| or L^4|D|/8' 
     233            ! 
     234            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90 
     235            ! 
     236            ! allocate locate gradient arrays used in ldf_dyn. Allocated once here for efficiency. 
     237            ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) ,  esqf(jpi,jpj) , STAT=ierr ) 
     238            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 
     239            ! 
     240            ! Set local gridscale values 
     241            DO jj = 2, jpjm1 
     242               DO ji = fs_2, fs_jpim1 
     243                  esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2  
     244                  esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2  
     245               END DO 
     246            END DO 
     247            ! 
    210248         CASE DEFAULT 
    211249            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 
     
    232270      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity)  
    233271      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator ) 
    234       !!                BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian 
    235272      !! 
    236       !! ** action  :    ahmt, ahmf   update at each time step 
     273      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)   
     274      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator ) 
     275      !! 
     276      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian 
     277      !! ** action  :    ahmt, ahmf   updated at each time step 
    237278      !!---------------------------------------------------------------------- 
    238279      INTEGER, INTENT(in) ::   kt   ! time step index 
     
    240281      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    241282      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax   ! local scalar 
     283      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb              ! local scalar 
    242284      !!---------------------------------------------------------------------- 
    243285      ! 
     
    248290      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity ) 
    249291         ! 
    250          IF( ln_dynldf_lap   ) THEN          !  laplacian operator : |u| e /12 = |u/144| e 
     292         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e 
    251293            DO jk = 1, jpkm1 
    252294               DO jj = 2, jpjm1 
     
    280322         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. ) 
    281323         ! 
     324         ! 
     325      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky) 
     326         ! 
     327         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D| 
     328            ! 
     329            zcmsmag = (rn_csmc/rpi)**2                                              ! (C_smag/pi)^2 
     330            zstabf_lo  = rn_cfacmin * rn_cfacmin / ( 2_.wp * 4._wp * zcmsmag )      !  lower limit stability factor scaling 
     331            zstabf_up  = rn_cfacmax / ( 4._wp * zcmsmag * rdt )                     !  upper limit stability factor scaling 
     332            ! 
     333            DO jk = 1, jpkm1 
     334               ! 
     335               DO jj = 2, jpj 
     336                  DO ji = 2, jpi 
     337                     zdb = ( (  ub(ji,jj,jk) * r1_e2u(ji,jj) -  ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) )  & 
     338                          &                  * r1_e1t(ji,jj) * e2t(ji,jj)                           & 
     339                          & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  & 
     340                          &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk) 
     341                     dtensq(ji,jj) = zdb*zdb 
     342                  END DO 
     343               END DO 
     344               ! 
     345               DO jj = 1, jpjm1 
     346                  DO ji = 1, jpim1 
     347                     zdb = ( (  ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) -  ub(ji,jj,jk) * r1_e1u(ji,jj) )  & 
     348                          &                    * r1_e2f(ji,jj)   * e1f(ji,jj)                       & 
     349                          & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  & 
     350                          &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk) 
     351                     dshesq(ji,jj) = zdb*zdb 
     352                  END DO 
     353               END DO 
     354               ! 
     355               DO jj = 2, jpjm1 
     356                  DO ji = fs_2, fs_jpim1 
     357                     ! 
     358                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk) 
     359                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
     360                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
     361                                                     ! T-point value 
     362                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
     363                     ahmt(ji,jj,jk) = zdelta * sqrt(          dtensq(ji,jj)   +                        & 
     364                                     &               r1_4 * ( dshesq(ji,jj)   + dshesq(ji,jj-1)   +    & 
     365                                     &                        dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 
     366                     ahmt(ji,jj,jk) =   MAX( ahmt(ji,jj,jk),   & 
     367                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == cfacmin * |U|L/2 
     368                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == cfacmax * L^2/(4dt) 
     369                                                     ! F-point value 
     370                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
     371                     ahmf(ji,jj,jk) = zdelta * sqrt(          dshesq(ji,jj)   +                        & 
     372                                     &               r1_4 * ( dtensq(ji,jj)   + dtensq(ji,jj+1)   +    & 
     373                                     &                        dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 
     374                     ahmf(ji,jj,jk) =   MAX( ahmf(ji,jj,jk),   & 
     375                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == cfacmin * |U|L/2 
     376                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == cfacmax * L^2/(4dt) 
     377                     ! 
     378                  END DO 
     379               END DO 
     380            END DO 
     381            ! 
     382         ENDIF 
     383         ! 
     384         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 
     385                                       !                      = sqrt( A_lap_smag L^2/8 ) 
     386                                       ! stability limits already applied to laplacian values 
     387            ! 
     388            DO jk = 1, jpkm1 
     389               ! 
     390               DO jj = 2, jpjm1 
     391                  DO ji = fs_2, fs_jpim1 
     392                     ! 
     393                     ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 
     394                     ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 
     395                     ! 
     396                  END DO 
     397               END DO 
     398            END DO 
     399            ! 
     400         ENDIF 
     401         ! 
     402         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. ) 
     403         ! 
    282404      END SELECT 
    283405      ! 
Note: See TracChangeset for help on using the changeset viewer.