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

Ignore:
Timestamp:
2017-03-17T10:44:05+01:00 (7 years ago)
Author:
jcastill
Message:

Changes as in HZG wave forcing branch, but adapted to r6232

File:
1 edited

Legend:

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

    r6204 r7807  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
     26   USE sbcwave 
     27   !  
    2628   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2729   USE lib_mpp        ! MPP manager 
     
    4749   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points 
    4850 
     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 
    4955   !                              !! ** Namelist  namzdf_gls  ** 
     56   LOGICAL  ::   ln_crban      = .FALSE. !! WAVE2NEMO added ln_crban (=.TRUE. use Craig and Banner scheme) 
    5057   LOGICAL  ::   ln_length_lim     ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 
    5158   LOGICAL  ::   ln_sigpsi         ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 
     
    5966   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2) 
    6067   REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 
    61    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 
    6269   REAL(wp) ::   rn_hsro           ! Minimum surface roughness 
    63    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)  
    6471 
    6572   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
     
    9299   REAL(wp) ::   rtrans        =  0.1_wp 
    93100   REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters 
    94    REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        - 
    95    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                  !     -           -           -        - 
    96103   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        - 
    97104   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
     
    115122      !!                ***  FUNCTION zdf_gls_alloc  *** 
    116123      !!---------------------------------------------------------------------- 
    117       ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    118          &      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),rsbc_tke1(jpi,jpj),        & 
     126         &      rsbc_tke3(jpi,jpj), rsbc_psi1(jpi,jpj)   , STAT=zdf_gls_alloc ) 
    119127         ! 
    120128      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    161169      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    162170 
     171 
     172      ! variable initialization 
     173      IF( ln_crban ) 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 
     177      ENDIF 
     178 
    163179      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    164180         avt (:,:,:) = avt_k (:,:,:) 
     
    197213         zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10) 
    198214         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    199       ! 
     215      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)  
     216         WHERE( swh_wavepar == 0._wp ) ! surface roughness length according to Charnock formula when sign. wave height 0 
     217            zhsro = MAX(rn_charn / grav * ustars2, rn_hsro) 
     218         ELSEWHERE 
     219            zhsro = MAX(swh_wavepar, rn_hsro) 
     220         END WHERE 
    200221      END SELECT 
    201222 
     
    314335      ! 
    315336      CASE ( 0 )             ! Dirichlet case 
    316       ! First level 
    317       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    318       en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    319       z_elem_a(:,:,1) = en(:,:,1) 
    320       z_elem_c(:,:,1) = 0._wp 
    321       z_elem_b(:,:,1) = 1._wp 
    322       !  
    323       ! One level below 
    324       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    325           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    326       en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    327       z_elem_a(:,:,2) = 0._wp  
    328       z_elem_c(:,:,2) = 0._wp 
    329       z_elem_b(:,:,2) = 1._wp 
    330       ! 
    331       ! 
     337         IF( ln_crban ) 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         ELSE                 ! wave induced mixing case with default values 
     350            en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
     351            en(:,:,1) = MAX(en(:,:,1), 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) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) & 
     358                &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
     359            en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
     360            z_elem_a(:,:,2) = 0._wp  
     361            z_elem_c(:,:,2) = 0._wp 
     362            z_elem_b(:,:,2) = 1._wp 
     363            ! 
     364            ! 
     365         ENDIF 
    332366      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    333       ! 
    334       ! Dirichlet conditions at k=1 
    335       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    336       en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
    337       z_elem_a(:,:,1) = en(:,:,1) 
    338       z_elem_c(:,:,1) = 0._wp 
    339       z_elem_b(:,:,1) = 1._wp 
    340       ! 
    341       ! at k=2, set de/dz=Fw 
    342       !cbr 
    343       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    344       z_elem_a(:,:,2) = 0._wp 
    345       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
    346       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    347            &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
    348  
    349       en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
    350       ! 
    351       ! 
     367         IF( ln_crban ) 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) 
     378         ELSE                  ! Shear free case: d(e)/dz=Fw with default values 
     379            ! Dirichlet conditions at k=1 
     380            en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp) 
     381            en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
     382            z_elem_a(:,:,1) = en(:,:,1) 
     383            z_elem_c(:,:,1) = 0._wp 
     384            z_elem_b(:,:,1) = 1._wp 
     385            ! 
     386            ! at k=2, set de/dz=Fw 
     387            !cbr 
     388            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     389            z_elem_a(:,:,2) = 0._wp 
     390            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
     391            zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
     392                 &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
     393 
     394            en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
     395            ! 
     396            ! 
     397         ENDIF 
    352398      END SELECT 
    353399 
     
    536582      ! 
    537583      CASE ( 0 )             ! Dirichlet boundary conditions 
    538       ! 
    539       ! Surface value 
    540       zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
    541       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    542       z_elem_a(:,:,1) = psi(:,:,1) 
    543       z_elem_c(:,:,1) = 0._wp 
    544       z_elem_b(:,:,1) = 1._wp 
    545       ! 
    546       ! One level below 
    547       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
    548       zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
    549       psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    550       z_elem_a(:,:,2) = 0._wp 
    551       z_elem_c(:,:,2) = 0._wp 
    552       z_elem_b(:,:,2) = 1._wp 
    553       !  
    554       ! 
     584         IF( ln_crban ) THEN   ! Wave induced mixing case 
     585                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
     586                               ! balance between the production and the 
     587                               ! 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            ! 
     605         ELSE                  ! Wave induced mixing case with default values 
     606            ! Surface value 
     607            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
     608            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     609            z_elem_a(:,:,1) = psi(:,:,1) 
     610            z_elem_c(:,:,1) = 0._wp 
     611            z_elem_b(:,:,1) = 1._wp 
     612            ! 
     613            ! One level below 
     614            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 
     615            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
     616            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     617            z_elem_a(:,:,2) = 0._wp 
     618            z_elem_c(:,:,2) = 0._wp 
     619            z_elem_b(:,:,2) = 1._wp 
     620            !  
     621            ! 
     622         ENDIF 
    555623      CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz 
    556       ! 
    557       ! Surface value: Dirichlet 
    558       zdep(:,:)       = zhsro(:,:) * rl_sf 
    559       psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
    560       z_elem_a(:,:,1) = psi(:,:,1) 
    561       z_elem_c(:,:,1) = 0._wp 
    562       z_elem_b(:,:,1) = 1._wp 
    563       ! 
    564       ! Neumann condition at k=2 
    565       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    566       z_elem_a(:,:,2) = 0._wp 
    567       ! 
    568       ! Set psi vertical flux at the surface: 
    569       zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
    570       zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
    571       zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
    572       zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
    573              & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
    574       zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
    575       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
    576  
    577       !    
    578       ! 
     624         IF( ln_crban ) 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            ! 
     642         ELSE                 ! Wave induced mixing case with default values 
     643            ! Surface value: Dirichlet 
     644            zdep(:,:)       = zhsro(:,:) * rl_sf 
     645            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
     646            z_elem_a(:,:,1) = psi(:,:,1) 
     647            z_elem_c(:,:,1) = 0._wp 
     648            z_elem_b(:,:,1) = 1._wp 
     649            ! 
     650            ! Neumann condition at k=2 
     651            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
     652            z_elem_a(:,:,2) = 0._wp 
     653            ! 
     654            ! Set psi vertical flux at the surface: 
     655            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
     656            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
     657            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
     658                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
     659            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
     660                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
     661            zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
     662            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
     663            !    
     664            ! 
     665         ENDIF 
    579666      END SELECT 
    580667 
     
    866953      !! 
    867954      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    868          &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    869          &            rn_crban, rn_charn, rn_frac_hs,        & 
     955         &            rn_clim_galp, ln_crban, ln_sigpsi, rn_hsro,      & 
     956         &            rn_crban_default, rn_charn, rn_frac_hs,& 
    870957         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    871958         &            nn_stab_func, nn_clos 
     
    895982         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot 
    896983         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi 
    897          WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
     984         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default 
     985         WRITE(numout,*) '      Craig and Banner scheme                       ln_crban       = ', ln_crban 
    898986         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    899987         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     
    909997 
    910998      !                                !* Check of some namelist values 
    911       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    912       IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 
    913       IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 
    914       IF( nn_stab_func  < 0 .OR. nn_stab_func  > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 
    915       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) 
     999      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' )  
     1000      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' )  
     1001      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' )  
     1002!      IF( nn_z0_met == 3 .AND. .NOT.ln_sdw ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_sdw=T' )  
     1003      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' )  
     1004      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
    9161005 
    9171006      SELECT CASE ( nn_clos )          !* set the parameters for the chosen closure 
     
    10851174               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
    10861175 
    1087       IF ( rn_crban==0._wp ) THEN 
     1176      IF( .NOT. ln_crban .AND. rn_crban_default==0._wp ) THEN 
    10881177         rl_sf = vkarmn 
    10891178      ELSE 
     
    11241213      rc03  = rc02 * rc0 
    11251214      rc04  = rc03 * rc0 
    1126       rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking 
    1127       rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking  
    1128       zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
    1129       rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer  
     1215      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 
     1216      rsbc_tke2         = rdt * rn_crban_default / rl_sf 
     1217      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
     1218      rtrans            = 0.2_wp / zcr 
    11301219      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    11311220      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.