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 13249 for NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2020-07-04T10:22:08+02:00 (4 years ago)
Author:
clem
Message:

merge with Jerome's branch NEMO_4.03_IODRAG

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/ZDF/zdftke.F90

    r13006 r13249  
    8383   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1) 
    8484   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean 
    85    REAL(wp) ::      rn_eice   ! =0 ON below sea-ice, =4 OFF when ice fraction > 1/4    
    8685   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    8786   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     87   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)    
    8888 
    8989   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
     
    199199      REAL(wp), DIMENSION(:,:,:) , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points) 
    200200      ! 
    201       INTEGER ::   ji, jj, jk              ! dummy loop arguments 
     201      INTEGER ::   ji, jj, jk                  ! dummy loop arguments 
    202202      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars 
    203203      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3 
    204204      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient 
    205       REAL(wp) ::   zbbrau, zri                ! local scalars 
    206       REAL(wp) ::   zfact1, zfact2, zfact3     !   -         - 
    207       REAL(wp) ::   ztx2  , zty2  , zcof       !   -         - 
    208       REAL(wp) ::   ztau  , zdif               !   -         - 
    209       REAL(wp) ::   zus   , zwlc  , zind       !   -         - 
    210       REAL(wp) ::   zzd_up, zzd_lw             !   -         - 
     205      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars 
     206      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      - 
     207      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      - 
     208      REAL(wp) ::   ztau  , zdif               !   -      - 
     209      REAL(wp) ::   zus   , zwlc  , zind       !   -      - 
     210      REAL(wp) ::   zzd_up, zzd_lw             !   -      - 
    211211      INTEGER , DIMENSION(jpi,jpj)     ::   imlc 
    212       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc, zfr_i 
     212      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3 
    213213      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw 
    214214      !!-------------------------------------------------------------------- 
    215215      ! 
    216       zbbrau = rn_ebb / rau0       ! Local constant initialisation 
    217       zfact1 = -.5_wp * rdt  
    218       zfact2 = 1.5_wp * rdt * rn_ediss 
    219       zfact3 = 0.5_wp       * rn_ediss 
     216      zbbrau  =  rn_ebb / rau0       ! Local constant initialisation 
     217      zbbirau =  3.75_wp / rau0 
     218      zfact1  = -0.5_wp * rdt  
     219      zfact2  =  1.5_wp * rdt * rn_ediss 
     220      zfact3  =  0.5_wp       * rn_ediss 
     221      ! 
     222      ! ice fraction considered for attenuation of langmuir & wave breaking 
     223      SELECT CASE ( nn_eice ) 
     224      CASE( 0 )   ;   zice_fra(:,:) = 0._wp 
     225      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp ) 
     226      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:) 
     227      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 
     228      END SELECT 
    220229      ! 
    221230      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    222231      !                     !  Surface/top/bottom boundary condition on tke 
    223232      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    224        
     233      ! 
    225234      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    226235         DO ji = fs_2, fs_jpim1   ! vector opt. 
    227             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     236            en(ji,jj,1) = MAX( rn_emin0, ( ( 1._wp - fr_i(ji,jj) ) * zbbrau + & 
     237               &                                     fr_i(ji,jj)   * zbbirau ) * taum(ji,jj) ) * tmask(ji,jj,1) 
    228238         END DO 
    229239      END DO 
     
    257267                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  & 
    258268                     &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  ) 
    259                   en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) & 
     269                  en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &      
    260270                     &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) 
    261271               END DO 
     
    295305            DO ji = fs_2, fs_jpim1   ! vector opt. 
    296306               zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    297                zfr_i(ji,jj) = ( 1._wp - 4._wp * fr_i(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    298                IF (zfr_i(ji,jj) < 0. ) zfr_i(ji,jj) = 0. 
     307               zus3(ji,jj) = ( 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok 
    299308            END DO 
    300309         END DO          
     
    302311            DO jj = 2, jpjm1 
    303312               DO ji = fs_2, fs_jpim1   ! vector opt. 
    304                   IF ( zfr_i(ji,jj) /= 0. ) THEN                
     313                  IF ( zus3(ji,jj) /= 0. ) THEN                
    305314                     ! vertical velocity due to LC    
    306315                     IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 
    307316                        !                                           ! vertical velocity due to LC 
    308                         zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )   ! warning: optimization: zus^3 is in zfr_i 
     317                        zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    309318                        !                                           ! TKE Langmuir circulation source term 
    310                         en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
     319                        en(ji,jj,jk) = en(ji,jj,jk) + rdt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) 
    311320                     ENDIF 
    312321                  ENDIF 
     
    408417       
    409418      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction) 
    410          DO jk = 2, jpkm1                       ! rn_eice =0 ON below sea-ice, =4 OFF when ice fraction > 0.25 
     419         DO jk = 2, jpkm1                       ! nn_eice=0 : ON below sea-ice ; nn_eice>0 : partly OFF 
    411420            DO jj = 2, jpjm1 
    412421               DO ji = fs_2, fs_jpim1   ! vector opt. 
    413422                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    414                      &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     423                     &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    415424               END DO 
    416425            END DO 
     
    421430               jk = nmln(ji,jj) 
    422431               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    423                   &                                 * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     432                  &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    424433            END DO 
    425434         END DO 
     
    434443                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    435444                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   & 
    436                      &                        * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
     445                     &                                 * ( 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    437446               END DO 
    438447            END DO 
     
    499508      zmxlm(:,:,:)  = rmxl_min     
    500509      zmxld(:,:,:)  = rmxl_min 
    501       ! 
    502      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
     510      !  
     511      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    503512         ! 
    504513         zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    505514#if ! defined key_si3 && ! defined key_cice 
    506          DO jj = 2, jpjm1 
     515         DO jj = 2, jpjm1                     ! No sea-ice 
    507516            DO ji = fs_2, fs_jpim1 
    508517               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1) 
     
    510519         END DO 
    511520#else 
     521 
    512522         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
    513523         ! 
     
    519529            END DO 
    520530            ! 
    521          CASE( 1 )                           ! scaling with constant sea-ice thickness 
     531         CASE( 1 )                      ! scaling with constant sea-ice thickness 
    522532            DO jj = 2, jpjm1 
    523533               DO ji = fs_2, fs_jpim1 
    524                   zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1) 
     534                  zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     535                     &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1) 
    525536               END DO 
    526537            END DO 
    527538            ! 
    528          CASE( 2 )                                 ! scaling with mean sea-ice thickness 
     539         CASE( 2 )                      ! scaling with mean sea-ice thickness 
    529540            DO jj = 2, jpjm1 
    530541               DO ji = fs_2, fs_jpim1 
    531542#if defined key_si3 
    532                   zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1) 
     543                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     544                     &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1) 
    533545#elif defined key_cice 
    534546                  zmaxice = MAXVAL( h_i(ji,jj,:) ) 
    535                   zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     547                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     548                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    536549#endif 
    537550               END DO 
    538551            END DO 
    539552            ! 
    540          CASE( 3 )                                 ! scaling with max sea-ice thickness 
     553         CASE( 3 )                      ! scaling with max sea-ice thickness 
    541554            DO jj = 2, jpjm1 
    542555               DO ji = fs_2, fs_jpim1 
    543556                  zmaxice = MAXVAL( h_i(ji,jj,:) ) 
    544                   zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1) 
     557                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + & 
     558                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1) 
    545559               END DO 
    546560            END DO 
     
    704718         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             & 
    705719         &                 nn_pdl  , ln_drg   , ln_lc    , rn_lc,      & 
    706          &                 nn_etau , nn_htau  , rn_efr   , rn_eice   
     720         &                 nn_etau , nn_htau  , rn_efr   , nn_eice   
    707721      !!---------------------------------------------------------------------- 
    708722      ! 
     
    737751            IF( nn_mxlice == 1 ) & 
    738752            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice 
     753            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice 
     754            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice' 
     755            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness' 
     756            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness' 
     757            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness' 
     758            CASE DEFAULT 
     759               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4') 
     760            END SELECT 
    739761         ENDIF 
    740762         WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg 
     
    744766         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau 
    745767         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr 
    746          WRITE(numout,*) '          below sea-ice:  =0 ON                      rn_eice   = ', rn_eice 
    747          WRITE(numout,*) '          =4 OFF when ice fraction > 1/4   ' 
     768         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice 
     769         SELECT CASE( nn_eice )  
     770         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking' 
     771         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )' 
     772         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)' 
     773         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )' 
     774         CASE DEFAULT 
     775            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 
     776         END SELECT       
    748777         IF( ln_drg ) THEN 
    749778            WRITE(numout,*) 
Note: See TracChangeset for help on using the changeset viewer.