Changeset 7809


Ignore:
Timestamp:
2017-03-17T15:21:06+01:00 (4 years ago)
Author:
jcastill
Message:

First implementation of the HZG wave focing/coupling branch - only ln_phioc and ln_tauoc in place. This is crashing amm7 runs with Baltic boundary conditions, but not without them.

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

Legend:

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

    r7797 r7809  
    269269   ln_sdw  = .false.       !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)    
    270270   ln_tauoc= .false.       !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)   
     271   ln_phioc= .false.       !  Activate wave to ocean energy (T => ln_wave=.true. & fill namsbc_wave)   
    271272   ln_stcor= .false.       !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)   
    272273   nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 
     
    10061007   rn_clim_galp  = 0.267   !  galperin limit 
    10071008   ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case 
    1008    rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux 
     1009   rn_crban_default = 100. !  Craig and Banner 1994 constant for wb tke flux 
    10091010   rn_charn      = 70000.  !  Charnock constant for wb induced roughness length 
    10101011   rn_hsro       =  0.02   !  Minimum surface roughness 
     
    13001301   sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
    13011302   sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
     1303   sn_phioc    =  'sdw_wave' ,        1          , 'wave_energy',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
    13021304!  
    13031305   cn_dir      = './'  !  root directory for the location of drag coefficient files 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r7792 r7809  
    6868   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
    6969   LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used  
     70   LOGICAL , PUBLIC ::   ln_phioc       !: true if wave energy to ocean is used  
    7071   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used  
    7172   ! 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r7797 r7809  
    11951195      !                                                      ! ========================= !   
    11961196      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1)  
     1197       
     1198      !                                                      ! ========================= !   
     1199      !                                                      !   Wave to ocean energy    !  
     1200      !                                                      ! ========================= !   
     1201      IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN 
     1202         rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1) 
     1203      ENDIF 
    11971204       
    11981205      !  Fields received by SAS when OASIS coupling 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r7792 r7809  
    9090         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
    9191         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx , nn_components, ln_cpl  ,   & 
    92          &             ln_wavcpl , nn_drag 
     92         &             ln_phioc  , ln_wavcpl , nn_drag 
    9393      INTEGER  ::   ios 
    9494      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm 
     
    136136         WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw  
    137137         WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc    = ', ln_tauoc  
     138         WRITE(numout,*) '                 wave to ocean energy - wave breaking    ln_phioc    = ', ln_phioc 
    138139         WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor    = ', ln_stcor  
    139140         WRITE(numout,*) '                 neutral drag coefficient                ln_cdgw     = ', ln_cdgw 
     
    224225      IF ( ln_wave ) THEN 
    225226         !Activated wave module but neither drag nor stokes drift activated  
    226          IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN   
    227              CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')  
     227         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc) )   THEN   
     228             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F, ln_phioc=F')  
    228229         ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN   
    229230             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
    230231         ENDIF 
    231232      ELSE 
    232          IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) &   
     233         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor .OR. ln_phioc ) &   
    233234            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &  
    234235            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &  
    235236            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &  
    236237            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      &  
    237             &                  'or Stokes-Coriolis term (ln_stcori=T)'  ) 
     238            &                  'or wave to ocean energy modification (ln_phioc=T) ',           &  
     239            &                  'or Stokes-Coriolis term (ln_stcor=T)'  ) 
    238240      ENDIF  
    239241      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r7797 r7809  
    5656   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao 
    5757   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean 
     58   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy  
    5859   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !: 
    5960   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !: 
     61   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing 
    6062   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !: 
    6163   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !: 
     
    222224      ENDIF 
    223225 
     226      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==! 
     227         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 
     229         WHERE( rn_crban <  10.0 ) rn_crban =  10.0 
     230         WHERE( rn_crban > 300.0 ) rn_crban = 300.0 
     231      ENDIF 
     232 
    224233      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==!  
    225234         ! 
     
    267276      CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    268277      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read 
    269       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  & 
     278      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  sn_phioc, & 
    270279                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read 
    271280      ! 
    272       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc 
     281      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc, sn_phioc 
    273282      !!--------------------------------------------------------------------- 
    274283      ! 
     
    285294         IF( .NOT. cpl_wdrag ) THEN 
    286295            ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    287             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     296            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' ) 
    288297            ! 
    289298                                   ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
     
    297306         IF( .NOT. cpl_tauoc ) THEN 
    298307            ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc 
    299             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     308            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' ) 
    300309            ! 
    301310                                    ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   ) 
     
    304313         ENDIF 
    305314         ALLOCATE( tauoc_wave(jpi,jpj) ) 
     315      ENDIF 
     316 
     317      IF( ln_phioc ) THEN 
     318         IF( .NOT. cpl_phioc ) THEN 
     319            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc 
     320            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' ) 
     321            ! 
     322                                    ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   ) 
     323            IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) ) 
     324            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' ) 
     325         ENDIF 
     326         ALLOCATE( rn_crban(jpi,jpj) ) 
    306327      ENDIF 
    307328 
     
    334355            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp 
    335356            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift 
    336             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' ) 
     357            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' ) 
    337358            ! 
    338359            DO ifpr= 1, jpfld 
     
    354375         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN 
    355376            ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum 
    356             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' ) 
     377            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' ) 
    357378                                   ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   ) 
    358379            IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) ) 
  • branches/UKMO/r6232_HZG_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r7481 r7809  
    2424   USE phycst         ! physical constants 
    2525   USE zdfmxl         ! mixed layer 
    26    USE sbcwave ,  ONLY: hsw   ! significant wave height  
     26   USE sbcwave, ONLY: hsw,rn_crban 
    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 
    5054 
    5155   !                              !! ** Namelist  namzdf_gls  ** 
     
    6165   REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2) 
    6266   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 
     67   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing 
    6468   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)  
     69   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1)  
    6670 
    6771   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
     
    9498   REAL(wp) ::   rtrans        =  0.1_wp 
    9599   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                  !     -           -           -        - 
     100   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        - 
     101   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        - 
    98102   REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        - 
    99103   REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        - 
     
    117121      !!                ***  FUNCTION zdf_gls_alloc  *** 
    118122      !!---------------------------------------------------------------------- 
    119       ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     & 
    120          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc ) 
     123      ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) ,     & 
     124         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     & 
     125         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     & 
     126         &      rsbc_psi1(jpi,jpj), STAT= zdf_gls_alloc ) 
    121127         ! 
    122128      IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc ) 
     
    163169      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    164170 
     171 
     172      ! variable initialization 
     173      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 
     177      ENDIF 
     178 
    165179      IF( kt /= nit000 ) THEN   ! restore before value to compute tke 
    166180         avt (:,:,:) = avt_k (:,:,:) 
     
    200214         zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 
    201215      CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)  
    202          zhsro(:,:) = hsw(:,:) 
     216         WHERE( hsw == 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(hsw, rn_hsro) 
     220         END WHERE 
    203221      END SELECT 
    204222 
     
    317335      ! 
    318336      CASE ( 0 )             ! Dirichlet case 
    319       ! First level 
    320       en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    321       en(:,:,1) = MAX(en(:,:,1), rn_emin)  
    322       z_elem_a(:,:,1) = en(:,:,1) 
    323       z_elem_c(:,:,1) = 0._wp 
    324       z_elem_b(:,:,1) = 1._wp 
    325       !  
    326       ! One level below 
    327       en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 
    328           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
    329       en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
    330       z_elem_a(:,:,2) = 0._wp  
    331       z_elem_c(:,:,2) = 0._wp 
    332       z_elem_b(:,:,2) = 1._wp 
    333       ! 
    334       ! 
     337         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         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 
    335366      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
    336       ! 
    337       ! Dirichlet conditions at k=1 
    338       en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 
    339       en(:,:,1)       = MAX(en(:,:,1), 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       ! at k=2, set de/dz=Fw 
    345       !cbr 
    346       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
    347       z_elem_a(:,:,2) = 0._wp 
    348       zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 
    349       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 
    350            &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
    351  
    352       en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
    353       ! 
    354       ! 
     367         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) 
     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 
    355398      END SELECT 
    356399 
     
    539582      ! 
    540583      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       ! 
     584         IF( ln_phioc ) 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 
    558623      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       ! 
     624         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            ! 
     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 
    582666      END SELECT 
    583667 
     
    870954      NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 
    871955         &            rn_clim_galp, ln_sigpsi, rn_hsro,      & 
    872          &            rn_crban, rn_charn, rn_frac_hs,        & 
     956         &            rn_crban_default, rn_charn, rn_frac_hs,& 
    873957         &            nn_bc_surf, nn_bc_bot, nn_z0_met,      & 
    874958         &            nn_stab_func, nn_clos 
     
    898982         WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot 
    899983         WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi 
    900          WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban 
     984         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default 
    901985         WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn 
    902986         WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met 
     
    915999      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' )  
    9161000      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' )  
     1001      IF( nn_z0_met == 3 .AND. .NOT.ln_sdw ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_sdw=T' )  
    9181002      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' )  
    9191003      IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' ) 
     
    10891173               &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 
    10901174 
    1091       IF ( rn_crban==0._wp ) THEN 
     1175      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN 
    10921176         rl_sf = vkarmn 
    10931177      ELSE 
     
    11281212      rc03  = rc02 * rc0 
    11291213      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  
     1214      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 
     1215      rsbc_tke2         = rdt * rn_crban_default / rl_sf 
     1216      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
     1217      rtrans            = 0.2_wp / zcr 
    11341218      rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness 
    11351219      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.