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 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90 – NEMO

Ignore:
Timestamp:
2017-02-06T10:25:03+01:00 (7 years ago)
Author:
timgraham
Message:

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90

    r6140 r7646  
    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_minfac       !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 
     48   REAL(wp), PUBLIC ::   rn_maxfac       !: 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_minfac, rn_maxfac 
    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_minfac    = ', rn_minfac 
     137         WRITE(numout,*) '         factor multiplier for theorectical lower upper for ' 
     138         WRITE(numout,*) '           Smagorinsky eddy visc (def. 1.0)   rn_maxfac    = ', rn_maxfac 
    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 arrays used in ldf_dyn.  
     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_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )        ! lower limit stability factor scaling 
     331            zstabf_up  = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )              ! upper limit stability factor scaling 
     332            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead  
     333            !                                                                       ! of |U|L^3/16 in blp case 
     334            DO jk = 1, jpkm1 
     335               ! 
     336               DO jj = 2, jpj 
     337                  DO ji = 2, jpi 
     338                     zdb = ( (  ub(ji,jj,jk) * r1_e2u(ji,jj) -  ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) )  & 
     339                          &                  * r1_e1t(ji,jj) * e2t(ji,jj)                           & 
     340                          & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  & 
     341                          &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk) 
     342                     dtensq(ji,jj) = zdb*zdb 
     343                  END DO 
     344               END DO 
     345               ! 
     346               DO jj = 1, jpjm1 
     347                  DO ji = 1, jpim1 
     348                     zdb = ( (  ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) -  ub(ji,jj,jk) * r1_e1u(ji,jj) )  & 
     349                          &                    * r1_e2f(ji,jj)   * e1f(ji,jj)                       & 
     350                          & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  & 
     351                          &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk) 
     352                     dshesq(ji,jj) = zdb*zdb 
     353                  END DO 
     354               END DO 
     355               ! 
     356               DO jj = 2, jpjm1 
     357                  DO ji = fs_2, fs_jpim1 
     358                     ! 
     359                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk) 
     360                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
     361                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
     362                                                     ! T-point value 
     363                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
     364                     ahmt(ji,jj,jk) = zdelta * sqrt(          dtensq(ji,jj)   +                        & 
     365                                     &               r1_4 * ( dshesq(ji,jj)   + dshesq(ji,jj-1)   +    & 
     366                                     &                        dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 
     367                     ahmt(ji,jj,jk) =   MAX( ahmt(ji,jj,jk),   & 
     368                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
     369                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt) 
     370                                                     ! F-point value 
     371                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
     372                     ahmf(ji,jj,jk) = zdelta * sqrt(          dshesq(ji,jj)   +                        & 
     373                                     &               r1_4 * ( dtensq(ji,jj)   + dtensq(ji,jj+1)   +    & 
     374                                     &                        dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 
     375                     ahmf(ji,jj,jk) =   MAX( ahmf(ji,jj,jk),   & 
     376                                     &   SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
     377                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt) 
     378                     ! 
     379                  END DO 
     380               END DO 
     381            END DO 
     382            ! 
     383         ENDIF 
     384         ! 
     385         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 
     386                                       !                      = sqrt( A_lap_smag L^2/8 ) 
     387                                       ! stability limits already applied to laplacian values 
     388                                       ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 
     389            ! 
     390            DO jk = 1, jpkm1 
     391               ! 
     392               DO jj = 2, jpjm1 
     393                  DO ji = fs_2, fs_jpim1 
     394                     ! 
     395                     ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 
     396                     ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 
     397                     ! 
     398                  END DO 
     399               END DO 
     400            END DO 
     401            ! 
     402         ENDIF 
     403         ! 
     404         CALL lbc_lnk( ahmt, 'T', 1. )   ;   CALL lbc_lnk( ahmf, 'F', 1. ) 
     405         ! 
    282406      END SELECT 
    283407      ! 
Note: See TracChangeset for help on using the changeset viewer.