Ignore:
Timestamp:
2019-03-20T19:36:44+01:00 (20 months ago)
Author:
clem
Message:

solve ticket #2261

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/LDF/ldfdyn.F90

    r10425 r10784  
    6262 
    6363   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 
    64    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,: ::   dtensq       !: horizontal tension squared         (Smagorinsky only) 
    65    REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,: ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only) 
     64   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only) 
     65   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only) 
    6666   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2            
    6767 
     
    242242         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 
    243243         ! 
    244          ahmt(:,:,jpk) = 0._wp                     ! last level always 0  
    245          ahmf(:,:,jpk) = 0._wp 
     244         ahmt(:,:,:) = 0._wp                       ! init to 0 needed  
     245         ahmf(:,:,:) = 0._wp 
    246246         ! 
    247247         !                                         ! value of lap/blp eddy mixing coef. 
     
    310310            ! 
    311311            !                          ! allocate arrays used in ldf_dyn.  
    312             ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 
     312            ALLOCATE( dtensq(jpi,jpj,jpk) , dshesq(jpi,jpj,jpk) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 
    313313            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 
    314314            ! 
    315             DO jj = 2, jpjm1           ! Set local gridscale values 
    316                DO ji = fs_2, fs_jpim1 
    317                   esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2  
    318                   esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2  
     315            DO jj = 1, jpj             ! Set local gridscale values 
     316               DO ji = 1, jpi 
     317                  esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2  
     318                  esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2  
    319319               END DO 
    320320            END DO 
     
    359359      ! 
    360360      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    361       REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax   ! local scalar 
    362       REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb              ! local scalar 
     361      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zemax   ! local scalar (option 31) 
     362      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb     ! local scalar (option 32) 
    363363      !!---------------------------------------------------------------------- 
    364364      ! 
     
    373373               DO jj = 2, jpjm1 
    374374                  DO ji = fs_2, fs_jpim1 
     375                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
     376                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
     377                     zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 
     378                     ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk)      ! 288= 12*12 * 2 
     379                  END DO 
     380               END DO 
     381               DO jj = 1, jpjm1 
     382                  DO ji = 1, fs_jpim1 
    375383                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk) 
    376384                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
    377                      zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
    378                      zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 
    379                      zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 
    380                      ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk)      ! 288= 12*12 * 2 
    381                      ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk) 
     385                     zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 
     386                     ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk)      ! 288= 12*12 * 2 
    382387                  END DO 
    383388               END DO 
     
    387392               DO jj = 2, jpjm1 
    388393                  DO ji = fs_2, fs_jpim1 
     394                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
     395                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
     396                     zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 
     397                     ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax  ) * zemax * tmask(ji,jj,jk) 
     398                  END DO 
     399               END DO 
     400               DO jj = 1, jpjm1 
     401                  DO ji = 1, fs_jpim1 
    389402                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk) 
    390403                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
    391                      zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
    392                      zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 
    393                      zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 
    394                      ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax  ) * zetmax * tmask(ji,jj,jk) 
    395                      ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax  ) * zefmax * fmask(ji,jj,jk) 
     404                     zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 
     405                     ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax  ) * zemax * fmask(ji,jj,jk) 
    396406                  END DO 
    397407               END DO 
     
    406416         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D| 
    407417            ! 
    408             zcmsmag = (rn_csmc/rpi)**2                                              ! (C_smag/pi)^2 
    409             zstabf_lo  = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )        ! lower limit stability factor scaling 
    410             zstabf_up  = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )              ! upper limit stability factor scaling 
     418            zcmsmag   = (rn_csmc/rpi)**2                                            ! (C_smag/pi)^2 
     419            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )         ! lower limit stability factor scaling 
     420            zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt )               ! upper limit stability factor scaling 
    411421            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead  
    412422            !                                                                       ! of |U|L^3/16 in blp case 
    413423            DO jk = 1, jpkm1 
    414424               ! 
    415                DO jj = 2, jpj 
    416                   DO ji = 2, jpi 
    417                      zdb = ( (  ub(ji,jj,jk) * r1_e2u(ji,jj) -  ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) )  & 
    418                           &                  * r1_e1t(ji,jj) * e2t(ji,jj)                           & 
    419                           & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) -  vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) )  & 
    420                           &                  * r1_e2t(ji,jj) * e1t(ji,jj)    ) * tmask(ji,jj,jk) 
    421                      dtensq(ji,jj) = zdb * zdb 
     425               DO jj = 2, jpjm1 
     426                  DO ji = 2, jpim1 
     427                     zdb =   ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj)  & 
     428                        &  - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj) 
     429                     dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk) 
    422430                  END DO 
    423431               END DO 
     
    425433               DO jj = 1, jpjm1 
    426434                  DO ji = 1, jpim1 
    427                      zdb = ( (  ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) -  ub(ji,jj,jk) * r1_e1u(ji,jj) )  & 
    428                           &                    * r1_e2f(ji,jj)   * e1f(ji,jj)                       & 
    429                           & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) -  vb(ji,jj,jk) * r1_e2v(ji,jj) )  & 
    430                           &                    * r1_e1f(ji,jj)   * e2f(ji,jj)  ) * fmask(ji,jj,jk) 
    431                      dshesq(ji,jj) = zdb * zdb 
     435                     zdb =   ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj)  & 
     436                        &  + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj) 
     437                     dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk) 
    432438                  END DO 
    433439               END DO 
    434440               ! 
    435                DO jj = 2, jpjm1 
     441            END DO 
     442            ! 
     443            CALL lbc_lnk_multi( 'ldfdyn', dtensq, 'T', 1. )  ! lbc_lnk on dshesq not needed 
     444            ! 
     445            DO jk = 1, jpkm1 
     446              ! 
     447               DO jj = 2, jpjm1                                ! T-point value 
    436448                  DO ji = fs_2, fs_jpim1 
     449                     ! 
     450                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
     451                     zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
     452                     ! 
     453                     zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
     454                     ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         & 
     455                        &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  & 
     456                        &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 
     457                     ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
     458                     ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt) 
     459                     ! 
     460                  END DO 
     461               END DO 
     462               ! 
     463               DO jj = 1, jpjm1                                ! F-point value 
     464                  DO ji = 1, fs_jpim1 
    437465                     ! 
    438466                     zu2pv2_ij_p1 = ub(ji  ,jj+1,jk) * ub(ji  ,jj+1,jk) + vb(ji+1,jj  ,jk) * vb(ji+1,jj  ,jk) 
    439467                     zu2pv2_ij    = ub(ji  ,jj  ,jk) * ub(ji  ,jj  ,jk) + vb(ji  ,jj  ,jk) * vb(ji  ,jj  ,jk) 
    440                      zu2pv2_ij_m1 = ub(ji-1,jj  ,jk) * ub(ji-1,jj  ,jk) + vb(ji  ,jj-1,jk) * vb(ji  ,jj-1,jk) 
    441                                                      ! T-point value 
    442                      zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
    443                      ahmt(ji,jj,jk) = zdelta * sqrt(          dtensq(ji,jj)   +                        & 
    444                                      &               r1_4 * ( dshesq(ji,jj)   + dshesq(ji,jj-1)   +    & 
    445                                      &                        dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 
    446                      ahmt(ji,jj,jk) =   MAX( ahmt(ji,jj,jk),   & 
    447                                      &   SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
    448                      ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt) 
    449                                                      ! F-point value 
     468                     ! 
    450469                     zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2 
    451                      ahmf(ji,jj,jk) = zdelta * sqrt(          dshesq(ji,jj)   +                        & 
    452                                      &               r1_4 * ( dtensq(ji,jj)   + dtensq(ji,jj+1)   +    & 
    453                                      &                        dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 
    454                      ahmf(ji,jj,jk) =   MAX( ahmf(ji,jj,jk),   & 
    455                                      &   SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
    456                      ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up )                    ! Impose upper limit == maxfac  * L^2/(4*2dt) 
     470                     ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         & 
     471                        &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  & 
     472                        &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 
     473                     ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2 
     474                     ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt) 
    457475                     ! 
    458476                  END DO 
    459477               END DO 
     478               ! 
    460479            END DO 
    461480            ! 
     
    470489                  DO ji = fs_2, fs_jpim1 
    471490                     ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 
     491                  END DO 
     492               END DO 
     493               DO jj = 1, jpjm1 
     494                  DO ji = 1, fs_jpim1 
    472495                     ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 
    473496                  END DO 
Note: See TracChangeset for help on using the changeset viewer.