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 10478 for branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 – NEMO

Ignore:
Timestamp:
2019-01-09T12:11:55+01:00 (5 years ago)
Author:
jcastill
Message:

Merged branch UKMO/r6232_HZG_WAVE-coupling@9868

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r10473 r10478  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
    26    USE sbcwave ,  ONLY: hsw   ! significant wave height   
     26   USE sbcwave 
    2727   ! 
    2828   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     
    4848   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points 
    4949   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points 
     50 
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke1  
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke3  
     53   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_psi1  
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rtrans 
    5055 
    5156   !                              !! ** Namelist  namzdf_gls  ** 
     
    6166   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2) 
    6267   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 
    63    REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing 
     68   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing 
    6469   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
    65    REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1)  
     70   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1, 3)  
    6671 
    6772   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
     
    9297   REAL(wp) ::   rm7           =  0.0_wp 
    9398   REAL(wp) ::   rm8           =  0.318_wp 
    94    REAL(wp) ::   rtrans        =  0.1_wp 
     99   REAL(wp) ::   rtrans_default =  0.1_wp 
    95100   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
    96    REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        - 
    97    REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        - 
     101   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        -  
     102   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        - 
    98103   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        - 
    99104   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
     
    117122      !!                ***  FUNCTION zdf_gls_alloc  *** 
    118123      !!---------------------------------------------------------------------- 
    119       ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    120          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
     124      ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) ,     &  
     125         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     &  
     126         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     &  
     127         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc ) 
    121128         ! 
    122129      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    163170      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    164171 
     172      ! variable initialization  
     173      IF( ln_phioc ) THEN  
     174         IF( nn_wmix==jp_janssen ) THEN  
     175            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking   
     176            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking  
     177            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking  
     178         ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     179            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf  
     180            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf  
     181            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp )  
     182         ENDIF  
     183      ENDIF  
     184 
    165185      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    166186         avt (:,:,:) = avt_k (:,:,:) 
     
    200220         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    201221      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)   
    202          zhsro(:,:) = hsw(:,:) 
     222         zhsro = MAX(rn_frac_hs * hsw, rn_hsro) 
    203223      END SELECT 
    204224 
     
    314334      ! --------------------------------- 
    315335      ! 
     336      IF( ln_phioc ) THEN  
     337         SELECT CASE ( nn_bc_surf )  
     338         !  
     339         CASE ( 0 )             ! Dirichlet case  
     340            IF( nn_wmix==jp_janssen ) THEN  
     341               ! First level  
     342               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )  
     343               z_elem_a(:,:,1) = en(:,:,1)  
     344               z_elem_c(:,:,1) = 0._wp  
     345               z_elem_b(:,:,1) = 1._wp  
     346  
     347               ! One level below  
     348               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin )  
     349               z_elem_a(:,:,2) = 0._wp  
     350               z_elem_c(:,:,2) = 0._wp  
     351               z_elem_b(:,:,2) = 1._wp  
     352            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     353               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)  
     354               en(:,:,1) = MAX(en(:,:,1), rn_emin)   
     355               z_elem_a(:,:,1) = en(:,:,1)  
     356               z_elem_c(:,:,1) = 0._wp  
     357               z_elem_b(:,:,1) = 1._wp  
     358               !   
     359               ! One level below  
     360               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) &  
     361                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)  
     362               en(:,:,2) = MAX(en(:,:,2), rn_emin )  
     363               z_elem_a(:,:,2) = 0._wp   
     364               z_elem_c(:,:,2) = 0._wp  
     365               z_elem_b(:,:,2) = 1._wp  
     366               !  
     367            ENDIF  
     368         CASE ( 1 )             ! Neumann boundary condition on d(e)/dz  
     369            IF( nn_wmix==jp_janssen ) THEN  
     370               ! Dirichlet conditions at k=1  
     371               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin )  
     372               z_elem_a(:,:,1) = en(:,:,1)  
     373               z_elem_c(:,:,1) = 0._wp  
     374               z_elem_b(:,:,1) = 1._wp  
     375               ! at k=2, set de/dz=Fw  
     376               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     377               z_elem_a(:,:,2) = 0._wp  
     378               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf)  
     379               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     380            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     381               ! Dirichlet conditions at k=1  
     382               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp)  
     383               en(:,:,1)       = MAX(en(:,:,1), rn_emin)        
     384               z_elem_a(:,:,1) = en(:,:,1)  
     385               z_elem_c(:,:,1) = 0._wp  
     386               z_elem_b(:,:,1) = 1._wp  
     387               !  
     388               ! at k=2, set de/dz=Fw  
     389               !cbr  
     390               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     391               z_elem_a(:,:,2) = 0._wp  
     392               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) ))  
     393               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) &  
     394                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)  
     395 
     396               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2)  
     397               !  
     398            ENDIF  
     399         END SELECT  
     400      ELSE 
    316401      SELECT CASE ( nn_bc_surf ) 
    317402      ! 
    318403      CASE ( 0 )             ! Dirichlet case 
    319404      ! First level 
    320       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     405      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
    321406      en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    322407      z_elem_a(:,:,1) = en(:,:,1) 
     
    325410      !  
    326411      ! One level below 
    327       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
     412      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    328413          &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    329414      en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
     
    334419      ! 
    335420      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    336       ! 
    337421      ! Dirichlet conditions at k=1 
    338       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
     422      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
    339423      en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
    340424      z_elem_a(:,:,1) = en(:,:,1) 
     
    346430      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    347431      z_elem_a(:,:,2) = 0._wp 
    348       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
     432      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) )) 
    349433      zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    350434           &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     
    354438      ! 
    355439      END SELECT 
     440      ENDIF ! ln_phioc 
    356441 
    357442      ! Bottom boundary condition on tke 
     
    539624      ! 
    540625      CASE ( 0 )             ! Dirichlet boundary conditions 
    541       ! 
    542       ! Surface value 
    543       zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
    544       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    545       z_elem_a(:,:,1) = psi(:,:,1) 
    546       z_elem_c(:,:,1) = 0._wp 
    547       z_elem_b(:,:,1) = 1._wp 
    548       ! 
    549       ! One level below 
    550       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
    551       zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
    552       psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    553       z_elem_a(:,:,2) = 0._wp 
    554       z_elem_c(:,:,2) = 0._wp 
    555       z_elem_b(:,:,2) = 1._wp 
    556       !  
    557       ! 
     626         IF( ln_phioc ) THEN   ! Wave induced mixing case  
     627                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2  
     628                               ! balance between the production and the  
     629                               ! dissipation terms including the wave effect  
     630            IF( nn_wmix==jp_janssen ) THEN  
     631               ! Surface value  
     632               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     633               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     634               z_elem_a(:,:,1) = psi(:,:,1)        
     635               z_elem_c(:,:,1) = 0._wp  
     636               z_elem_b(:,:,1) = 1._wp  
     637               !  
     638               ! One level below  
     639               zex1 = (rmm*ra_sf+rnn)  
     640               zex2 = (rmm*ra_sf)  
     641               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2  
     642               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1)  
     643               z_elem_a(:,:,2) = 0._wp  
     644               z_elem_c(:,:,2) = 0._wp  
     645               z_elem_b(:,:,2) = 1._wp  
     646               !   
     647               !  
     648            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     649               ! Surface value  
     650               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     651               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     652               z_elem_a(:,:,1) = psi(:,:,1)  
     653               z_elem_c(:,:,1) = 0._wp  
     654               z_elem_b(:,:,1) = 1._wp  
     655               !  
     656               ! One level below  
     657               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) )))  
     658               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)  
     659               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     660               z_elem_a(:,:,2) = 0._wp  
     661               z_elem_c(:,:,2) = 0._wp  
     662               z_elem_b(:,:,2) = 1._wp  
     663               !   
     664               !  
     665            ENDIF  
     666         ELSE                  ! Wave induced mixing case with default values  
     667            ! Surface value  
     668            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic  
     669            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     670            z_elem_a(:,:,1) = psi(:,:,1)  
     671            z_elem_c(:,:,1) = 0._wp  
     672            z_elem_b(:,:,1) = 1._wp  
     673            !  
     674            ! One level below  
     675            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) )))  
     676            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)  
     677            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     678            z_elem_a(:,:,2) = 0._wp  
     679            z_elem_c(:,:,2) = 0._wp  
     680            z_elem_b(:,:,2) = 1._wp  
     681            !   
     682            !  
     683         ENDIF 
    558684      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    559       ! 
    560       ! Surface value: Dirichlet 
    561       zdep(:,:)       = zhsro(:,:) * rl_sf 
    562       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    563       z_elem_a(:,:,1) = psi(:,:,1) 
    564       z_elem_c(:,:,1) = 0._wp 
    565       z_elem_b(:,:,1) = 1._wp 
    566       ! 
    567       ! Neumann condition at k=2 
    568       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    569       z_elem_a(:,:,2) = 0._wp 
    570       ! 
    571       ! Set psi vertical flux at the surface: 
    572       zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    573       zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    574       zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    575       zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
    576              & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
    577       zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
    578       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    579  
    580       !    
    581       ! 
     685         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields  
     686            IF( nn_wmix==jp_janssen ) THEN  
     687               !  
     688               zdep(:,:) = rl_sf * zhsro(:,:)  
     689               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     690               z_elem_a(:,:,1) = psi(:,:,1)  
     691               z_elem_c(:,:,1) = 0._wp  
     692               z_elem_b(:,:,1) = 1._wp  
     693               !  
     694               ! Neumann condition at k=2  
     695               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     696               z_elem_a(:,:,2) = 0._wp  
     697               !  
     698               ! Set psi vertical flux at the surface:  
     699               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf)  
     700               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) &  
     701                  &                   * en(:,:,1)**rmm * zdep  
     702               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     703               !  
     704            ELSE IF( nn_wmix==jp_craigbanner ) THEN  
     705               ! Surface value: Dirichlet  
     706               zdep(:,:)       = zhsro(:,:) * rl_sf  
     707               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     708               z_elem_a(:,:,1) = psi(:,:,1)  
     709               z_elem_c(:,:,1) = 0._wp  
     710               z_elem_b(:,:,1) = 1._wp  
     711               !  
     712               ! Neumann condition at k=2  
     713               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     714               z_elem_a(:,:,2) = 0._wp  
     715               !  
     716               ! Set psi vertical flux at the surface:  
     717               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope  
     718               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)  
     719               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * &  
     720                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)  
     721               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &  
     722                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)  
     723               zflxs(:,:) = zdep(:,:) * zflxs(:,:)  
     724               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     725               !     
     726            ENDIF  
     727         ELSE                 ! Wave induced mixing case with default values  
     728            ! Surface value: Dirichlet  
     729            zdep(:,:)       = zhsro(:,:) * rl_sf  
     730            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)  
     731            z_elem_a(:,:,1) = psi(:,:,1)  
     732            z_elem_c(:,:,1) = 0._wp  
     733            z_elem_b(:,:,1) = 1._wp  
     734            !  
     735            ! Neumann condition at k=2  
     736            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b  
     737            z_elem_a(:,:,2) = 0._wp  
     738            !  
     739            ! Set psi vertical flux at the surface:  
     740            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope  
     741            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)  
     742            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * &  
     743                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp)  
     744            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &  
     745                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)  
     746            zflxs(:,:) = zdep(:,:) * zflxs(:,:)  
     747            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)  
     748            !     
     749            !  
     750         ENDIF 
    582751      END SELECT 
    583752 
     
    8701039      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    8711040         &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    872          &            rn_crban, rn_charn, rn_frac_hs,        & 
     1041         &            rn_crban_default, rn_charn, rn_frac_hs,& 
    8731042         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    8741043         &            nn_stab_func, nn_clos 
     
    8981067         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot 
    8991068         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi 
    900          WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
     1069         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default 
    9011070         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    9021071         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     
    9151084      IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_bc_surf is 0 or 1' )   
    9161085      IF( nn_z0_met < 0 .OR. nn_z0_met > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_z0_met is 0, 1, 2 or 3' )   
    917       IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )   
     1086      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' )   
     1087      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN  
     1088         CALL ctl_warn('W A R N I N G:  ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3')   
     1089         nn_z0_met = 3   
     1090      ENDIF 
    9181091      IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_stab_func is 0, 1, 2 and 3' )   
    9191092      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
     
    10891262               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
    10901263 
    1091       IF ( rn_crban==0._wp ) THEN 
     1264      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 
    10921265         rl_sf = vkarmn 
    10931266      ELSE 
     
    11281301      rc03  = rc02 * rc0 
    11291302      rc04  = rc03 * rc0 
    1130       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1131       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
    1132       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1133       rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1303      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf  
     1304      rsbc_tke2         = rdt * rn_crban_default / rl_sf  
     1305      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp )  
     1306      rtrans_default    = 0.2_wp / zcr 
    11341307      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    11351308      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.