Changeset 7853


Ignore:
Timestamp:
2017-03-30T16:22:25+02:00 (3 years ago)
Author:
jcastill
Message:

Add the original Craig and Banner vertical mixing scheme in case of wave coupling - further checks on wave forcing fields in case the ocean and wave land/sea masks are not the same

Location:
branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref

    r7809 r7853  
    10171017   nn_stab_func  =     2   !  stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB) 
    10181018   nn_clos       =     1   !  predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen) 
     1019   nn_wmix       =     0   !  type of wave breaking mixing 
     1020           !                             ! = 0 standard formulation (original NEMO formulation)  
     1021           !                             ! = 1 Craig and Banner formulation  
    10191022/ 
    10201023!----------------------------------------------------------------------- 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7809 r7853  
    226226      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==! 
    227227         CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing 
    228          rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rw/ra)sbc_wa 
     228         rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density 
    229229         WHERE( rn_crban <  10.0 ) rn_crban =  10.0 
    230230         WHERE( rn_crban > 300.0 ) rn_crban = 300.0 
     
    235235         IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled 
    236236            CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing 
    237             IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
    238             IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
    239             IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
    240             IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     237            IF( jp_hsw > 0 ) THEN 
     238               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height 
     239               WHERE( hsw > 100.0 ) hsw = 0.0 
     240               WHERE( hsw <   0.0 ) hsw = 0.0 
     241            ENDIF 
     242            IF( jp_wmp > 0 ) THEN 
     243               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period 
     244               WHERE( wmp > 100.0 ) wmp = 0.0 
     245               WHERE( wmp <   0.0 ) wmp = 0.0 
     246            ENDIF 
     247            IF( jp_usd > 0 ) THEN 
     248               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point 
     249               WHERE( ut0sd < -100.0 ) ut0sd = 1.0 
     250               WHERE( ut0sd >  100.0 ) ut0sd = 1.0 
     251            ENDIF 
     252            IF( jp_vsd > 0 ) THEN 
     253               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point 
     254               WHERE( vt0sd < -100.0 ) vt0sd = 1.0 
     255               WHERE( vt0sd >  100.0 ) vt0sd = 1.0 
     256            ENDIF 
    241257         ENDIF 
    242258         ! 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7809 r7853  
    5252   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_tke3 
    5353   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rsbc_psi1 
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rtrans 
    5455 
    5556   !                              !! ** Namelist  namzdf_gls  ** 
     
    6162   INTEGER  ::   nn_stab_func      ! stability functions G88, KC or Canuto (=0/1/2) 
    6263   INTEGER  ::   nn_clos           ! closure 0/1/2/3 MY82/k-eps/k-w/gen 
     64   INTEGER  ::   nn_wmix           ! type of wave breaking mixing 
     65   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0   ! Craig and Banner formulation (original NEMO formulation - 
     66                                                        !    direct conversion of mechanical to turbulent energy) 
     67   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1   ! Janssen formulation - no assumption on direct energy conversion  
    6368   REAL(wp) ::   rn_clim_galp      ! Holt 2008 value for k-eps: 0.267 
    6469   REAL(wp) ::   rn_epsmin         ! minimum value of dissipation (m2/s3) 
     
    96101   REAL(wp) ::   rm7           =  0.0_wp 
    97102   REAL(wp) ::   rm8           =  0.318_wp 
    98    REAL(wp) ::   rtrans        =  0.1_wp 
     103   REAL(wp) ::   rtrans_default =  0.1_wp 
    99104   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
    100105   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        - 
     
    124129         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     & 
    125130         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     & 
    126          &      rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc ) 
     131         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc ) 
    127132         ! 
    128133      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    172177      ! variable initialization 
    173178      IF( ln_phioc ) THEN 
    174          rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking  
    175          rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking 
    176          rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking 
     179         IF( nn_wmix==jp_janssen ) THEN 
     180            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking  
     181            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking 
     182            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking 
     183         ELSE 
     184            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf 
     185            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf 
     186            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
     187         ENDIF 
    177188      ENDIF 
    178189 
     
    336347      CASE ( 0 )             ! Dirichlet case 
    337348         IF( ln_phioc ) THEN  ! wave induced mixing case with forced/coupled fields 
    338             ! First level 
    339             en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
    340             z_elem_a(:,:,1) = en(:,:,1) 
    341             z_elem_c(:,:,1) = 0._wp 
    342             z_elem_b(:,:,1) = 1._wp 
    343  
    344             ! One level below 
    345             en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
    346             z_elem_a(:,:,2) = 0._wp 
    347             z_elem_c(:,:,2) = 0._wp 
    348             z_elem_b(:,:,2) = 1._wp 
     349            IF( nn_wmix==jp_janssen ) THEN 
     350               ! First level 
     351               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
     352               z_elem_a(:,:,1) = en(:,:,1) 
     353               z_elem_c(:,:,1) = 0._wp 
     354               z_elem_b(:,:,1) = 1._wp 
     355 
     356               ! One level below 
     357               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
     358               z_elem_a(:,:,2) = 0._wp 
     359               z_elem_c(:,:,2) = 0._wp 
     360               z_elem_b(:,:,2) = 1._wp 
     361            ELSE 
     362               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 
     363               en(:,:,1) = MAX(en(:,:,1), rn_emin)  
     364               z_elem_a(:,:,1) = en(:,:,1) 
     365               z_elem_c(:,:,1) = 0._wp 
     366               z_elem_b(:,:,1) = 1._wp 
     367               !  
     368               ! One level below 
     369               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) & 
     370                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
     371               en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
     372               z_elem_a(:,:,2) = 0._wp  
     373               z_elem_c(:,:,2) = 0._wp 
     374               z_elem_b(:,:,2) = 1._wp 
     375               ! 
     376            ENDIF 
    349377         ELSE                 ! wave induced mixing case with default values 
    350378            en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
     
    366394      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    367395         IF( ln_phioc ) THEN   ! Shear free case: d(e)/dz=Fw with forced/coupled fields 
    368             ! Dirichlet conditions at k=1 
    369             en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
    370             z_elem_a(:,:,1) = en(:,:,1) 
    371             z_elem_c(:,:,1) = 0._wp 
    372             z_elem_b(:,:,1) = 1._wp 
    373             ! at k=2, set de/dz=Fw 
    374             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    375             z_elem_a(:,:,2) = 0._wp 
    376             zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 
    377             en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     396            IF( nn_wmix==jp_janssen ) THEN 
     397               ! Dirichlet conditions at k=1 
     398               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
     399               z_elem_a(:,:,1) = en(:,:,1) 
     400               z_elem_c(:,:,1) = 0._wp 
     401               z_elem_b(:,:,1) = 1._wp 
     402               ! at k=2, set de/dz=Fw 
     403               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     404               z_elem_a(:,:,2) = 0._wp 
     405               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 
     406               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     407            ELSE 
     408               ! Dirichlet conditions at k=1 
     409               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 
     410               en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
     411               z_elem_a(:,:,1) = en(:,:,1) 
     412               z_elem_c(:,:,1) = 0._wp 
     413               z_elem_b(:,:,1) = 1._wp 
     414               ! 
     415               ! at k=2, set de/dz=Fw 
     416               !cbr 
     417               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     418               z_elem_a(:,:,2) = 0._wp 
     419               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) )) 
     420               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) & 
     421                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     422 
     423               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
     424               ! 
     425            ENDIF 
    378426         ELSE                  ! Shear free case: d(e)/dz=Fw with default values 
    379427            ! Dirichlet conditions at k=1 
     
    388436            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    389437            z_elem_a(:,:,2) = 0._wp 
    390             zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
     438            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) )) 
    391439            zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    392440                 &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     
    586634                               ! balance between the production and the 
    587635                               ! dissipation terms including the wave effect 
    588             ! Surface value 
    589             zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
    590             psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    591             z_elem_a(:,:,1) = psi(:,:,1)       
    592             z_elem_c(:,:,1) = 0._wp 
    593             z_elem_b(:,:,1) = 1._wp 
    594             ! 
    595             ! One level below 
    596             zex1 = (rmm*ra_sf+rnn) 
    597             zex2 = (rmm*ra_sf) 
    598             zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 
    599             psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
    600             z_elem_a(:,:,2) = 0._wp 
    601             z_elem_c(:,:,2) = 0._wp 
    602             z_elem_b(:,:,2) = 1._wp 
    603             !  
    604             ! 
     636            IF( nn_wmix==jp_janssen ) THEN 
     637               ! Surface value 
     638               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
     639               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     640               z_elem_a(:,:,1) = psi(:,:,1)       
     641               z_elem_c(:,:,1) = 0._wp 
     642               z_elem_b(:,:,1) = 1._wp 
     643               ! 
     644               ! One level below 
     645               zex1 = (rmm*ra_sf+rnn) 
     646               zex2 = (rmm*ra_sf) 
     647               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 
     648               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
     649               z_elem_a(:,:,2) = 0._wp 
     650               z_elem_c(:,:,2) = 0._wp 
     651               z_elem_b(:,:,2) = 1._wp 
     652               !  
     653               ! 
     654            ELSE 
     655               ! Surface value 
     656               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
     657               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     658               z_elem_a(:,:,1) = psi(:,:,1) 
     659               z_elem_c(:,:,1) = 0._wp 
     660               z_elem_b(:,:,1) = 1._wp 
     661               ! 
     662               ! One level below 
     663               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) ))) 
     664               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
     665               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     666               z_elem_a(:,:,2) = 0._wp 
     667               z_elem_c(:,:,2) = 0._wp 
     668               z_elem_b(:,:,2) = 1._wp 
     669               !  
     670               ! 
     671            ENDIF 
    605672         ELSE                  ! Wave induced mixing case with default values 
    606673            ! Surface value 
     
    612679            ! 
    613680            ! One level below 
    614             zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
     681            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) ))) 
    615682            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
    616683            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     
    623690      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    624691         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields 
    625             ! 
    626             zdep(:,:) = rl_sf * zhsro(:,:) 
    627             psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    628             z_elem_a(:,:,1) = psi(:,:,1) 
    629             z_elem_c(:,:,1) = 0._wp 
    630             z_elem_b(:,:,1) = 1._wp 
    631             ! 
    632             ! Neumann condition at k=2 
    633             z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    634             z_elem_a(:,:,2) = 0._wp 
    635             ! 
    636             ! Set psi vertical flux at the surface: 
    637             zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 
    638             zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
    639                &                   * en(:,:,1)**rmm * zdep 
    640             psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    641             ! 
     692            IF( nn_wmix==jp_janssen ) THEN 
     693               ! 
     694               zdep(:,:) = rl_sf * zhsro(:,:) 
     695               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     696               z_elem_a(:,:,1) = psi(:,:,1) 
     697               z_elem_c(:,:,1) = 0._wp 
     698               z_elem_b(:,:,1) = 1._wp 
     699               ! 
     700               ! Neumann condition at k=2 
     701               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     702               z_elem_a(:,:,2) = 0._wp 
     703               ! 
     704               ! Set psi vertical flux at the surface: 
     705               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 
     706               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
     707                  &                   * en(:,:,1)**rmm * zdep 
     708               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     709               ! 
     710            ELSE 
     711               ! Surface value: Dirichlet 
     712               zdep(:,:)       = zhsro(:,:) * rl_sf 
     713               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     714               z_elem_a(:,:,1) = psi(:,:,1) 
     715               z_elem_c(:,:,1) = 0._wp 
     716               z_elem_b(:,:,1) = 1._wp 
     717               ! 
     718               ! Neumann condition at k=2 
     719               z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     720               z_elem_a(:,:,2) = 0._wp 
     721               ! 
     722               ! Set psi vertical flux at the surface: 
     723               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
     724               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
     725               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
     726                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     727               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
     728                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
     729               zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
     730               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     731               !    
     732            ENDIF 
    642733         ELSE                 ! Wave induced mixing case with default values 
    643734            ! Surface value: Dirichlet 
     
    653744            ! 
    654745            ! Set psi vertical flux at the surface: 
    655             zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
     746            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    656747            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    657748            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
     
    9561047         &            rn_crban_default, rn_charn, rn_frac_hs,& 
    9571048         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    958          &            nn_stab_func, nn_clos 
     1049         &            nn_stab_func, nn_clos, nn_wmix 
    9591050      !!---------------------------------------------------------- 
    9601051      ! 
     
    12151306      rsbc_tke2         = rdt * rn_crban_default / rl_sf 
    12161307      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1217       rtrans            = 0.2_wp / zcr 
     1308      rtrans_default    = 0.2_wp / zcr 
    12181309      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    12191310      rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness  
Note: See TracChangeset for help on using the changeset viewer.