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 10922 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90 – NEMO

Ignore:
Timestamp:
2019-05-02T17:10:39+02:00 (5 years ago)
Author:
acc
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert IOM, LDF, OBS and SBC directories and compatibility changes elsewhere that these changes enforce. Changes pass SETTE and compare with original trunk results. Outstanding issues (currently with work-arounds) in DIU/step_diu.F90 and fld_bdy_interp within SBC/fldread.F90; proper soltions pending

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcisf.F90

    r10536 r10922  
    7575CONTAINS 
    7676  
    77   SUBROUTINE sbc_isf( kt ) 
     77  SUBROUTINE sbc_isf( kt, Kmm ) 
    7878      !!--------------------------------------------------------------------- 
    7979      !!                  ***  ROUTINE sbc_isf  *** 
    … …  
    8989      !!---------------------------------------------------------------------- 
    9090      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     91      INTEGER, INTENT(in) ::   Kmm  ! ocean time level indices 
    9192      ! 
    9293      INTEGER ::   ji, jj, jk   ! loop index 
    … …  
    102103         CASE ( 1 )    ! realistic ice shelf formulation 
    103104            ! compute T/S/U/V for the top boundary layer 
    104             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    105             CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T') 
    106             CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U') 
    107             CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V') 
     105            CALL sbc_isf_tbl(ts(:,:,:,jp_tem,Kmm),ttbl(:,:),'T',Kmm) 
     106            CALL sbc_isf_tbl(ts(:,:,:,jp_sal,Kmm),stbl(:,:),'T',Kmm) 
     107            CALL sbc_isf_tbl(uu(:,:,:,Kmm)       ,utbl(:,:),'U',Kmm) 
     108            CALL sbc_isf_tbl(vv(:,:,:,Kmm)       ,vtbl(:,:),'V',Kmm) 
    108109            ! iom print 
    109110            CALL iom_put('ttbl',ttbl(:,:)) 
    … …  
    113114            ! compute fwf and heat flux 
    114115            ! compute fwf and heat flux 
    115             IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt) 
     116            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt, Kmm) 
    116117            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfusisf  ! heat        flux 
    117118            ENDIF 
    … …  
    119120         CASE ( 2 )    ! Beckmann and Goosse parametrisation  
    120121            stbl(:,:)   = soce 
    121             CALL sbc_isf_bg03(kt) 
     122            CALL sbc_isf_bg03(kt, Kmm) 
    122123            ! 
    123124         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    … …  
    179180                  ikb = misfkb(ji,jj) 
    180181                  DO jk = ikt, ikb - 1 
    181                      zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
    182                      zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
    183                      zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk) 
     182                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm) 
     183                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm) 
     184                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm) 
    184185                  END DO 
    185186                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   &  
    186                      &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
     187                     &                                                                   * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm) 
    187188                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   &  
    188                      &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
     189                     &                                                                   * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm) 
    189190                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   &   
    190                      &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk) 
     191                     &                                                                   * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm) 
    191192               END DO 
    192193            END DO 
    … …  
    251252 
    252253 
    253   SUBROUTINE sbc_isf_init 
     254  SUBROUTINE sbc_isf_init( Kmm ) 
    254255      !!--------------------------------------------------------------------- 
    255256      !!                  ***  ROUTINE sbc_isf_init  *** 
    … …  
    263264      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
    264265      !!---------------------------------------------------------------------- 
     266      INTEGER, INTENT(in) ::   Kmm  ! ocean time level indices 
    265267      INTEGER               :: ji, jj, jk           ! loop index 
    266268      INTEGER               :: ik                   ! current level index 
    … …  
    355357                ik = 2 
    356358!!gm potential bug: use gdepw_0 not _n 
    357                 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
     359                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw(ji,jj,ik,Kmm) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
    358360                misfkt(ji,jj) = ik-1 
    359361            END DO 
    … …  
    386388            ikb = misfkt(ji,jj) 
    387389            ! thickness of boundary layer at least the top level thickness 
    388             rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 
     390            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t(ji,jj,ikt,Kmm)) 
    389391 
    390392            ! determine the deepest level influenced by the boundary layer 
    391393            DO jk = ikt+1, mbkt(ji,jj) 
    392                IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk 
    393             END DO 
    394             rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 
     394               IF( (SUM(e3t(ji,jj,ikt:jk-1,Kmm)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk 
     395            END DO 
     396            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t(ji,jj,ikt:ikb,Kmm))) ! limit the tbl to water thickness. 
    395397            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl 
    396398            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    397399 
    398             zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
    399             ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
     400            zhk           = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
     401            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t(ji,jj,ikb,Kmm)  ! proportion of bottom cell influenced by boundary layer 
    400402         END DO 
    401403      END DO 
    … …  
    411413 
    412414 
    413   SUBROUTINE sbc_isf_bg03(kt) 
     415  SUBROUTINE sbc_isf_bg03( kt, Kmm ) 
    414416      !!--------------------------------------------------------------------- 
    415417      !!                  ***  ROUTINE sbc_isf_bg03  *** 
    … …  
    426428      !!---------------------------------------------------------------------- 
    427429      INTEGER, INTENT ( in ) :: kt 
     430      INTEGER, INTENT ( in ) :: Kmm  ! ocean time level indices 
    428431      ! 
    429432      INTEGER  :: ji, jj, jk ! dummy loop index 
    … …  
    444447               DO jk = misfkt(ji,jj),misfkb(ji,jj) 
    445448                  ! Calculate freezing temperature 
    446                   zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04 
     449                  zpress = grav*rau0*gdept(ji,jj,ik,Kmm)*1.e-04 
    447450                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress)  
    448                   zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
     451                  zt_sum = zt_sum + (ts(ji,jj,jk,jp_tem,Kmm)-zt_frz) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)  ! sum temp 
    449452               END DO 
    450453               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
    … …  
    466469 
    467470 
    468   SUBROUTINE sbc_isf_cav( kt ) 
     471  SUBROUTINE sbc_isf_cav( kt, Kmm ) 
    469472      !!--------------------------------------------------------------------- 
    470473      !!                     ***  ROUTINE sbc_isf_cav  *** 
    … …  
    480483      !!--------------------------------------------------------------------- 
    481484      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     485      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    482486      ! 
    483487      INTEGER  ::   ji, jj     ! dummy loop indices 
    … …  
    520524 
    521525            ! compute gammat every where (2d) 
    522             CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     526            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, Kmm) 
    523527             
    524528            ! compute upward heat flux zhtflx and upward water flux zwflx 
    … …  
    536540         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
    537541            ! compute gammat every where (2d) 
    538             CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     542            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, Kmm) 
    539543 
    540544            ! compute upward heat flux zhtflx and upward water flux zwflx 
    … …  
    600604 
    601605 
    602    SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 
     606   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf, Kmm ) 
    603607      !!---------------------------------------------------------------------- 
    604608      !! ** Purpose    : compute the coefficient echange for heat flux 
    … …  
    611615      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !  
    612616      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !  
     617      INTEGER                 , INTENT(in   ) ::   Kmm  ! ocean time level indices 
    613618      ! 
    614619      INTEGER  :: ji, jj                     ! loop index 
    … …  
    679684!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation 
    680685!!gm moreover, use Max(rn2,0) to take care of static instabilities.... 
    681                   zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1) 
     686                  zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm) 
    682687                  !                                            ! shear of horizontal velocity 
    683                   zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
    684                      &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
    685                   zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  & 
    686                      &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
     688                  zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  & 
     689                     &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  ) 
     690                  zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  & 
     691                     &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  ) 
    687692                  !                                            ! richardson number (minimum value set to zero) 
    688693                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
    … …  
    691696                  zts(jp_tem) = ttbl(ji,jj) 
    692697                  zts(jp_sal) = stbl(ji,jj) 
    693                   zdep        = gdepw_n(ji,jj,ikt) 
     698                  zdep        = gdepw(ji,jj,ikt,Kmm) 
    694699                  ! 
    695700                  CALL eos_rab( zts, zdep, zab ) 
    … …  
    700705                  !! compute Monin Obukov Length 
    701706                  ! Maximum boundary layer depth 
    702                   zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp 
     707                  zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp 
    703708                  ! Compute Monin obukhov length scale at the surface and Ekman depth: 
    704709                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 
    … …  
    727732 
    728733 
    729    SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 
     734   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin, Kmm ) 
    730735      !!---------------------------------------------------------------------- 
    731736      !!                  ***  SUBROUTINE sbc_isf_tbl  *** 
    … …  
    737742      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout 
    738743      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out 
     744      INTEGER                   , INTENT(in   ) :: Kmm  ! ocean time level indices 
    739745      ! 
    740746      INTEGER ::   ji, jj, jk                ! loop index 
    … …  
    753759               ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
    754760               ! thickness of boundary layer at least the top level thickness 
    755                zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) ) 
     761               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u(ji,jj,ikt,Kmm) ) 
    756762 
    757763               ! determine the deepest level influenced by the boundary layer 
    758764               DO jk = ikt+1, mbku(ji,jj) 
    759                   IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
    760                END DO 
    761                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     765                  IF ( (SUM(e3u(ji,jj,ikt:jk-1,Kmm)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     766               END DO 
     767               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u(ji,jj,ikt:ikb,Kmm)))  ! limit the tbl to water thickness. 
    762768 
    763769               ! level fully include in the ice shelf boundary layer 
    764770               DO jk = ikt, ikb - 1 
    765                   ze3 = e3u_n(ji,jj,jk) 
     771                  ze3 = e3u(ji,jj,jk,Kmm) 
    766772                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
    767773               END DO 
    768774 
    769775               ! level partially include in ice shelf boundary layer  
    770                zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     776               zhk = SUM( e3u(ji, jj, ikt:ikb - 1,Kmm)) / zhisf_tbl(ji,jj) 
    771777               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    772778            END DO 
    … …  
    785791               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
    786792               ! thickness of boundary layer at least the top level thickness 
    787                zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt)) 
     793               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v(ji,jj,ikt,Kmm)) 
    788794 
    789795               ! determine the deepest level influenced by the boundary layer 
    790796               DO jk = ikt+1, mbkv(ji,jj) 
    791                   IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
    792                END DO 
    793                zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     797                  IF ( (SUM(e3v(ji,jj,ikt:jk-1,Kmm)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     798               END DO 
     799               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v(ji,jj,ikt:ikb,Kmm)))  ! limit the tbl to water thickness. 
    794800 
    795801               ! level fully include in the ice shelf boundary layer 
    796802               DO jk = ikt, ikb - 1 
    797                   ze3 = e3v_n(ji,jj,jk) 
     803                  ze3 = e3v(ji,jj,jk,Kmm) 
    798804                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
    799805               END DO 
    800806 
    801807               ! level partially include in ice shelf boundary layer  
    802                zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     808               zhk = SUM( e3v(ji, jj, ikt:ikb - 1,Kmm)) / zhisf_tbl(ji,jj) 
    803809               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    804810            END DO 
    … …  
    820826               ! level fully include in the ice shelf boundary layer 
    821827               DO jk = ikt, ikb - 1 
    822                   ze3 = e3t_n(ji,jj,jk) 
     828                  ze3 = e3t(ji,jj,jk,Kmm) 
    823829                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    824830               END DO 
    825831 
    826832               ! level partially include in ice shelf boundary layer  
    827                zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
     833               zhk = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj) 
    828834               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
    829835            END DO 
    … …  
    837843       
    838844 
    839    SUBROUTINE sbc_isf_div( phdivn ) 
     845   SUBROUTINE sbc_isf_div( phdivn, Kmm ) 
    840846      !!---------------------------------------------------------------------- 
    841847      !!                  ***  SUBROUTINE sbc_isf_div  *** 
    … …  
    850856      !!---------------------------------------------------------------------- 
    851857      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence 
     858      INTEGER                   , INTENT( in    ) ::   Kmm      ! ocean time level indices 
    852859      !  
    853860      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    865872               ikb = misfkt(ji,jj) 
    866873               ! thickness of boundary layer at least the top level thickness 
    867                rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt)) 
     874               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t(ji,jj,ikt,Kmm)) 
    868875 
    869876               ! determine the deepest level influenced by the boundary layer 
    870877               DO jk = ikt, mbkt(ji,jj) 
    871                   IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    872                END DO 
    873                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     878                  IF ( (SUM(e3t(ji,jj,ikt:jk-1,Kmm)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     879               END DO 
     880               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t(ji,jj,ikt:ikb,Kmm)))  ! limit the tbl to water thickness. 
    874881               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl 
    875882               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    876883 
    877                zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
    878                ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
     884               zhk           = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
     885               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t(ji,jj,ikb,Kmm)  ! proportion of bottom cell influenced by boundary layer 
    879886            END DO 
    880887         END DO 
Note: See TracChangeset for help on using the changeset viewer.