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 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 – NEMO

Ignore:
Timestamp:
2015-04-13T15:08:59+02:00 (9 years ago)
Author:
davestorkey
Message:

Merge in changes from trunk up to 5021.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r4624 r5208  
    241241      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242242         DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     243            IF (mikt(ji,jj) .GT. 1) THEN 
     244               en(ji,jj,mikt(ji,jj))=rn_emin 
     245            ELSE 
     246               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     247            END IF 
    244248         END DO 
    245249      END DO 
    … …  
    291295         END DO 
    292296         !                               ! finite LC depth 
    293 # if defined key_vectopt_loop 
    294          DO jj = 1, 1 
    295             DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    296 # else 
    297297         DO jj = 1, jpj  
    298298            DO ji = 1, jpi 
    299 # endif 
    300299               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) 
    301300            END DO 
    302301         END DO 
    303302         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
    304 !CDIR NOVERRCHK 
    305303         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    306 !CDIR NOVERRCHK 
    307304            DO jj = 2, jpjm1 
    308 !CDIR NOVERRCHK 
    309305               DO ji = fs_2, fs_jpim1   ! vector opt. 
    310306                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    … …  
    342338      END DO 
    343339      ! 
    344       DO jk = 2, jpkm1           !* Matrix and right hand side in en 
    345          DO jj = 2, jpjm1 
    346             DO ji = fs_2, fs_jpim1   ! vector opt. 
     340      DO jj = 2, jpjm1 
     341         DO ji = fs_2, fs_jpim1   ! vector opt. 
     342            DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en 
    347343               zcof   = zfact1 * tmask(ji,jj,jk) 
    348344               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
    … …  
    360356               !                                   ! right hand side in en 
    361357               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    362                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk) 
    363             END DO 
    364          END DO 
    365       END DO 
    366       !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    367       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    368          DO jj = 2, jpjm1 
    369             DO ji = fs_2, fs_jpim1    ! vector opt. 
     358                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     359                  &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     360            END DO 
     361            !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
     362            DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    370363               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    371364            END DO 
    372          END DO 
    373       END DO 
    374       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    375          DO ji = fs_2, fs_jpim1    ! vector opt. 
    376             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    377          END DO 
    378       END DO 
    379       DO jk = 3, jpkm1 
    380          DO jj = 2, jpjm1 
    381             DO ji = fs_2, fs_jpim1    ! vector opt. 
     365            ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     366            zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke 
     367            ! 
     368            DO jk = mikt(ji,jj)+2, jpkm1 
    382369               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    383370            END DO 
    384          END DO 
    385       END DO 
    386       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    387          DO ji = fs_2, fs_jpim1    ! vector opt. 
     371            ! 
     372            ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    388373            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    389          END DO 
    390       END DO 
    391       DO jk = jpk-2, 2, -1 
    392          DO jj = 2, jpjm1 
    393             DO ji = fs_2, fs_jpim1    ! vector opt. 
     374            ! 
     375            DO jk = jpk-2, mikt(ji,jj)+1, -1 
    394376               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    395377            END DO 
    396          END DO 
    397       END DO 
    398       DO jk = 2, jpkm1                             ! set the minimum value of tke 
    399          DO jj = 2, jpjm1 
    400             DO ji = fs_2, fs_jpim1   ! vector opt. 
     378            ! 
     379            DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke 
    401380               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
    402381            END DO 
    … …  
    412391               DO ji = fs_2, fs_jpim1   ! vector opt. 
    413392                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    414                      &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     393                     &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
    415394               END DO 
    416395            END DO 
    … …  
    421400               jk = nmln(ji,jj) 
    422401               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    423                   &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     402                  &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
    424403            END DO 
    425404         END DO 
    … …  
    433412                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    434413                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    435                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress  
     414                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
    436415                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    437416                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    438417                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    439                      &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 
     418                     &                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 
    440419               END DO 
    441420            END DO 
    … …  
    506485      ! 
    507486      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    508          zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    509          zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  ) 
    510       ELSE                          ! surface set to the minimum value 
    511          zmxlm(:,:,1) = rn_mxl0 
     487         DO jj = 2, jpjm1 
     488            DO ji = fs_2, fs_jpim1 
     489               IF (mikt(ji,jj) .GT. 1) THEN 
     490                  zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 
     491               ELSE 
     492                  zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     493                  zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 
     494               END IF 
     495            END DO 
     496         END DO 
     497      ELSE  
     498         DO jj = 2, jpjm1 
     499            DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value 
     500               zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 
     501            END DO 
     502         END DO 
    512503      ENDIF 
    513504      zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value 
    514505      ! 
    515506!CDIR NOVERRCHK 
    516       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    517 !CDIR NOVERRCHK 
    518          DO jj = 2, jpjm1 
    519 !CDIR NOVERRCHK 
    520             DO ji = fs_2, fs_jpim1   ! vector opt. 
     507      DO jj = 2, jpjm1 
     508!CDIR NOVERRCHK 
     509         DO ji = fs_2, fs_jpim1   ! vector opt. 
     510            !CDIR NOVERRCHK 
     511            DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    521512               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    522513               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    523514            END DO 
     515            zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value  
    524516         END DO 
    525517      END DO 
    … …  
    533525      ! 
    534526      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    535          DO jk = 2, jpkm1 
    536             DO jj = 2, jpjm1 
    537                DO ji = fs_2, fs_jpim1   ! vector opt. 
    538                   zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   & 
     527         DO jj = 2, jpjm1 
     528            DO ji = fs_2, fs_jpim1   ! vector opt. 
     529               DO jk = mikt(ji,jj)+1, jpkm1 
     530                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    539531                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    540532                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    545537         ! 
    546538      CASE ( 1 )           ! bounded by the vertical scale factor 
    547          DO jk = 2, jpkm1 
    548             DO jj = 2, jpjm1 
    549                DO ji = fs_2, fs_jpim1   ! vector opt. 
     539         DO jj = 2, jpjm1 
     540            DO ji = fs_2, fs_jpim1   ! vector opt. 
     541               DO jk = mikt(ji,jj)+1, jpkm1 
    550542                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    551543                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    556548         ! 
    557549      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    558          DO jk = 2, jpkm1         ! from the surface to the bottom : 
    559             DO jj = 2, jpjm1 
    560                DO ji = fs_2, fs_jpim1   ! vector opt. 
     550         DO jj = 2, jpjm1 
     551            DO ji = fs_2, fs_jpim1   ! vector opt. 
     552               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : 
    561553                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    562554               END DO 
    563             END DO 
    564          END DO 
    565          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
    566             DO jj = 2, jpjm1 
    567                DO ji = fs_2, fs_jpim1   ! vector opt. 
     555               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : 
    568556                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    569557                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    574562         ! 
    575563      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    576          DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
    577             DO jj = 2, jpjm1 
    578                DO ji = fs_2, fs_jpim1   ! vector opt. 
     564         DO jj = 2, jpjm1 
     565            DO ji = fs_2, fs_jpim1   ! vector opt. 
     566               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup 
    579567                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    580568               END DO 
    581             END DO 
    582          END DO 
    583          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
    584             DO jj = 2, jpjm1 
    585                DO ji = fs_2, fs_jpim1   ! vector opt. 
     569               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown 
    586570                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    587571               END DO 
    … …  
    628612      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    629613      ! 
    630       DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points 
     614      DO jj = 2, jpjm1 
     615         DO ji = fs_2, fs_jpim1   ! vector opt. 
     616            DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points 
     617               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
     618            END DO 
     619            DO jk = mikv(ji,jj)+1, jpkm1 
     620               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     621            END DO 
     622         END DO 
     623      END DO 
     624      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     625      ! 
     626      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    631627         DO jj = 2, jpjm1 
    632628            DO ji = fs_2, fs_jpim1   ! vector opt. 
    633                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    634                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
    635             END DO 
    636          END DO 
    637       END DO 
    638       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    639       ! 
    640       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    641          DO jk = 2, jpkm1 
    642             DO jj = 2, jpjm1 
    643                DO ji = fs_2, fs_jpim1   ! vector opt. 
     629               DO jk = mikt(ji,jj)+1, jpkm1 
    644630                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    645631                  !                                          ! shear 
    … …  
    655641                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    656642# if defined key_c1d 
    657                   e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number 
    658                   e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
     643                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
     644                  e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri 
    659645# endif 
    660646              END DO 
    … …  
    740726      ! 
    741727      !                               !* Check of some namelist values 
    742       IF( nn_mxl  < 0  .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
    743       IF( nn_pdl  < 0  .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
    744       IF( nn_htau < 0  .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
    745 #if ! key_coupled 
    746       IF( nn_etau == 3 )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    747 #endif 
     728      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
     729      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
     730      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     731      IF( nn_etau == 3 .AND. .NOT. lk_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 
    748732 
    749733      IF( ln_mxl0 ) THEN 
Note: See TracChangeset for help on using the changeset viewer.