Changeset 7807


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

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

Location:
branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r6204 r7807  
    297297              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    298298              ln_tra_dmp(ib_bdy)=.false. 
    299            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
    300               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    301299           ELSE 
    302300              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6204 r7807  
    767767 
    768768      IF( lk_vvl ) THEN 
    769          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content 
    770          CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content 
     769         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T  )   ! heat content 
     770         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T  )   ! salt content 
    771771         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content 
    772772         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r5407 r7807  
    6262   !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step 
    6363   !                                             !:  = 2 annual global mean of e-p-r set to zero 
    64    LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model 
    65    LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model 
    66    LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
     64  ! LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model 
     65   LOGICAL , PUBLIC ::   ln_cdgw = .FALSE.        !: true if neutral drag coefficient from wave model 
     66  ! LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
     67         ! WAVE2NEMO 11.12.2016 
     68      LOGICAL , PUBLIC ::   ln_wavetke  = .FALSE.   !: true if wave parametersare read 
     69      LOGICAL , PUBLIC ::   ln_stcor    = .FALSE.   !: true if Stokes-Coriolisforcing is included  
     70      LOGICAL , PUBLIC ::   ln_tauoc    = .FALSE.   !: true if wave-modifiedwater-side stress is used 
     71      ! END  WAVE2NEMO 
     72    
     73 
    6774   ! 
    6875   LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs 
     
    125132#endif 
    126133   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
    127  
     134   ! START WAVE2NEMO 26.11.2016 
     135   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tauoc_wavepar !:normalized stress into the ocean from wave model  CCC 
     136   REAL(wp), PUBLIC,ALLOCATABLE,SAVE,DIMENSION (:,:)       :: wspd_wavepar  !CCCmoved to sbc_oce    
     137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_wave !:wave model drag coefficient [N/m2] 
     138  ! END  
     139!WAVE2NEMO (14.12.2016) rn_crban moved to here, it needs to be a more global variable 
     140   REAL(wp),PUBLIC,ALLOCATABLE,SAVE, DIMENSION (:,:)  ::  rn_crban 
    128141   !!---------------------------------------------------------------------- 
    129142   !!                     Sea Surface Mean fields 
     
    176189      ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) ) 
    177190#endif 
    178          ! 
     191      ! Start WAVE2NEMO 11.12.2016 
     192      ALLOCATE( cdn_wave(jpi,jpj) ) 
     193      ALLOCATE( wspd_wavepar(jpi,jpj) ) 
     194      ! End 
     195   ! 
     196          ALLOCATE( rn_crban(jpi,jpj) ) !WAVE2NEMO (13.12.2016) 
     197 
     198! 
    179199      sbc_oce_alloc = MAXVAL( ierr ) 
    180200      IF( lk_mpp            )   CALL mpp_sum ( sbc_oce_alloc ) 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5582 r7807  
    1717   !!            3.7  !  2014-06  (L. Brodeau) simplification and optimization of CORE bulk 
    1818   !!---------------------------------------------------------------------- 
    19  
     19   !!            3.6     2016-06  (V.Alari) Added wave dependent drag 
     20   !coefficient and ocean side stress 
    2021   !!---------------------------------------------------------------------- 
    2122   !!   sbc_blk_core    : bulk formulation as ocean surface boundary condition (forced mode, CORE bulk formulea) 
     
    4142   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    4243   USE prtctl          ! Print control 
    43    USE sbcwave, ONLY   :  cdn_wave ! wave module 
     44!   USE sbcwave, ONLY   :  cdn_wave ! wave module (V.Alari 27.06.2016) We get 
     45!   info from sbc_oce 
    4446   USE sbc_ice         ! Surface boundary condition: ice fields 
    4547   USE lib_fortran     ! to use key_nosignedzero 
     
    183185         ! 
    184186         lhftau = ln_taudif                        ! do we use HF tau information? 
     187 
     188 
     189 
    185190         jfld = jpfld - COUNT( (/.NOT. lhftau/) ) 
    186191         ! 
     
    192197         END DO 
    193198         !                                         ! fill sf with slf_i and control print 
     199 
    194200         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 
    195201         ! 
     
    258264      REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height 
    259265      !!--------------------------------------------------------------------- 
     266      REAL(wp) ::   ztheta               ! local variable, wind direction,V.Alari 27.06.2016 
     267 
    260268      ! 
    261269      IF( nn_timing == 1 )  CALL timing_start('blk_oce_core') 
     
    285293      END DO 
    286294#endif 
     295        ! Start 27.06.2016 (V.Alari) 
     296           ! Wind vector at T points 
     297      IF ( ln_cdgw ) THEN 
     298         ! Use neutral 10-m wind speed from wave model if available 
     299         DO jj = 2, jpjm1 
     300            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     301               ! Local wind direction from non-neutral 10-m wind vector 
     302               ! (relative to grid and no current) 
     303               ! since we do not have the directional information from the wave 
     304               ! model 
     305               ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1),sf(jp_wndj)%fnow(ji,jj,1)) 
     306               ! Wind vector magnitude from 10-m neutral wind speed from wave 
     307               ! model 
     308               zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta) 
     309               zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta) 
     310               ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0 
     311               zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj  )+pu(ji,jj))) 
     312               zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji  ,jj-1)+pv(ji,jj))) 
     313            END DO 
     314         END DO 
     315      ELSE ! If no neutral wind from wave model 
    287316      DO jj = 2, jpjm1 
    288317         DO ji = fs_2, fs_jpim1   ! vect. opt. 
     
    291320         END DO 
    292321      END DO 
     322    ENDIF ! End (V.Alari) 
     323 
    293324      CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. ) 
    294325      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) 
     
    320351     
    321352      ! ... tau module, i and j component 
     353! Start 27.06.2016 (V.Alari) 
     354     IF ( ln_cdgw ) THEN 
     355         DO jj = 1, jpj 
     356            DO ji = 1, jpi 
     357               !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj) 
     358               zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj) 
     359               !taum(ji,jj) = tauoc_wavepar(ji,jj) 
     360               !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj) 
     361               taum  (ji,jj) = zztmp * wndm  (ji,jj) 
     362               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     363               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     364            ENDDO 
     365         ENDDO 
     366      ELSE 
    322367      DO jj = 1, jpj 
    323368         DO ji = 1, jpi 
     
    328373         END DO 
    329374      END DO 
    330  
     375     ENDIF ! End, (V.Alari) 
    331376      ! ... add the HF tau contribution to the wind stress module? 
    332377      IF( lhftau ) THEN 
     
    349394      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
    350395 
    351      
     396        ! START: 27.06.2015 (V. Alari). Added tauoc 
     397 
     398         IF (ln_tauoc) THEN 
     399         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 
     400         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 
     401         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 
     402         ENDIF 
     403 
     404        ! END  (V.Alari) 
    352405      !  Turbulent fluxes over ocean 
    353406      ! ----------------------------- 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_mfs.F90

    r5215 r7807  
    2424   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2525   USE prtctl          ! Print control 
    26    USE sbcwave,ONLY : cdn_wave !wave module 
     26!   USE sbcwave,ONLY : cdn_wave !wave module 
    2727 
    2828   IMPLICIT NONE 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r4990 r7807  
    2222   USE lib_mpp         ! distribued memory computing library 
    2323   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     24   USE sbcwave         ! Wave parameters 
    2425 
    2526   IMPLICIT NONE 
     
    152153               zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
    153154               zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
     155               IF (ln_cdgw) THEN 
     156                  zcdrag = 1.5e-3 
     157                  IF (cdn_wave(ji,jj)>=0.) zcdrag = cdn_wave(ji,jj) 
     158                  zcoef = 1./(zrhoa*zcdrag) 
     159               ENDIF 
    154160               taum(ji,jj) = zmod 
    155161               wndm(ji,jj) = SQRT( zmod * zcoef ) 
     
    158164         taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1) 
    159165         CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
     166 
     167         IF (ln_tauoc) THEN 
     168         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 
     169         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 
     170         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 
     171         ENDIF 
    160172 
    161173         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked) 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r5628 r7807  
    8585      !!---------------------------------------------------------------------- 
    8686      INTEGER ::   icpt   ! local integer 
    87       !! 
     87      !! V.Alari 26.06.2016 Added logicals for waves, deleted ln_wave ln_sdw, which came 
     88      !with UKMO 
    8889      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   & 
    8990         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   & 
    90          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
    91          &             nn_lsm    , nn_limflx , nn_components, ln_cpl 
     91         &             ln_ssr    , nn_isf    , nn_fwb,                                        & 
     92         &             nn_lsm    , nn_limflx , nn_components, ln_cpl,ln_stcor, ln_wavetke,    & 
     93         &             ln_tauoc, ln_cdgw 
    9294      INTEGER  ::   ios 
    9395      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm 
     
    120122          nn_ice      =   0 
    121123      ENDIF 
    122  
     124     ! V.Alari added wave logicalc (26.06.2016) 
    123125      IF(lwp) THEN               ! Control print 
    124126         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)' 
     
    146148         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea 
    147149         WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm 
    148       ENDIF 
     150         WRITE(numout,*) '         Stokes-Coriolis forcing from wavemodel ln_stcor       = ', ln_stcor 
     151         WRITE(numout,*) '              TKE wave breaking BC from wave model ln_wavetke        = ', ln_wavetke 
     152         WRITE(numout,*) '              Wave-modified stress from wave model ln_tauoc          = ', ln_tauoc 
     153         WRITE(numout,*) '              Drag coefficient from wave model ln_cdgw              = ', ln_cdgw 
     154  
     155 
     156     ENDIF 
    149157 
    150158      ! LIM3 Multi-category heat flux formulation 
     
    179187 
    180188      !                          ! Checks: 
    181       IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf  
     189      IF( nn_isf .EQ. 0 ) THEN                      ! no specific treatment in vicinity of ice shelf  
    182190         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 
    183191         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp 
     
    214222      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   & 
    215223         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' ) 
    216        
    217       IF ( ln_wave ) THEN 
     224    
     225     ! V.Alari comment out UKMO wave stuff    
     226    !  IF ( ln_wave ) THEN 
    218227      !Activated wave module but neither drag nor stokes drift activated 
    219          IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN 
    220             CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
     228     !    IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN 
     229      !      CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' ) 
    221230      !drag coefficient read from wave model definable only with mfs bulk formulae and core  
    222          ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
    223              CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
    224          ENDIF 
    225       ELSE 
    226       IF ( ln_cdgw .OR. ln_sdw  )                                         &  
    227          &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     & 
    228          & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 
    229       ENDIF  
     231       !  ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN        
     232        !     CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core') 
     233       !  ENDIF 
     234     ! ELSE 
     235     ! IF ( ln_cdgw .OR. ln_sdw  )                                         &  
     236      !   &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     & 
     237       !  & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 
     238     ! ENDIF  
     239      ! End V.Alari 
    230240      !                          ! Choice of the Surface Boudary Condition (set nsbc) 
    231241      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl 
     
    352362      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 
    353363      !                                                        ! averaged over nf_sbc time-step 
    354  
    355       IF (ln_wave) CALL sbc_wave( kt ) 
     364            ! START V.Alari 26.06.2016 
     365     IF (ln_cdgw) CALL sbc_wave( kt ) 
     366      ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for 
     367      ! TKE are needed 
     368      ! Also read it if drag coefficient to ensure that we use the neutral 10m 
     369      ! from WAM. 
     370      IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc) THEN 
     371            CALL sbc_wavepar( kt ) 
     372      ENDIF 
     373 
     374       ! END 
     375 
     376!      IF (ln_wave) CALL sbc_wave( kt ) 
    356377                                                   !==  sbc formulation  ==! 
    357378                                                             
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

    r5215 r7807  
    44   !! Wave module  
    55   !!====================================================================== 
    6    !! History :  3.3.1  !   2011-09  (Adani M)  Original code: Drag Coefficient  
    7    !!         :  3.4    !   2012-10  (Adani M)                 Stokes Drift  
     6   !! History :  3.3.1  !   2011-09  (Adani M)  Original code 
    87   !!---------------------------------------------------------------------- 
    98   USE iom             ! I/O manager library 
    109   USE in_out_manager  ! I/O manager 
    1110   USE lib_mpp         ! distribued memory computing library 
    12    USE fldread        ! read input fields 
    13    USE oce 
    14    USE sbc_oce        ! Surface boundary condition: ocean fields 
    15    USE domvvl 
    16  
    17     
     11   USE fldread         ! read input fields 
     12   USE sbc_oce         ! Surface boundary condition: ocean fields 
     13   USE phycst          ! physical constants 
     14   ! USE sbcblk_core, ONLY: Cd_n10    
    1815   !!---------------------------------------------------------------------- 
    1916   !!   sbc_wave       : read drag coefficient from wave model in netcdf files  
     
    2320   PRIVATE 
    2421 
    25    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
    26     
    27    INTEGER , PARAMETER ::   jpfld  = 3           ! maximum number of files to read for srokes drift 
    28    INTEGER , PARAMETER ::   jp_usd = 1           ! index of stokes drift  (i-component) (m/s)    at T-point 
    29    INTEGER , PARAMETER ::   jp_vsd = 2           ! index of stokes drift  (j-component) (m/s)    at T-point 
    30    INTEGER , PARAMETER ::   jp_wn  = 3           ! index of wave number                 (1/m)    at T-point 
    31    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient 
    32    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd    ! structure of input fields (file informations, fields read) Stokes Drift 
    33    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave  
    34    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum  
    35    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d  
     22   PUBLIC   sbc_wavepar    ! routine called in sbc_blk_core or sbc_blk_mfs 
     23    PUBLIC   sbc_wave    ! routine called in sbc_blk_core or sbc_blk_mfs 
     24 
     25   INTEGER , PARAMETER ::   jpfld_wavepar = 7    ! maximum number of fields to be read  
     26   INTEGER , PARAMETER ::   jp_ust = 1           ! index of Stokes velocity (east component) (m/s)  at T-point 
     27   INTEGER , PARAMETER ::   jp_vst = 2           ! index of Stokes velocity (north component) (m/s) at T-point 
     28   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave height (m) 
     29   INTEGER , PARAMETER ::   jp_mwp = 4           ! index of mean wave period (s) 
     30   INTEGER , PARAMETER ::   jp_phioc = 5         ! index of normalized energy flux into the ocean (non-dim) 
     31   INTEGER , PARAMETER ::   jp_tauoc = 6         ! index of normalized wave stress into the ocean (non-dim) 
     32   !INTEGER , PARAMETER ::   jp_cdww = 7          ! index of wave drag coeff (non-dim) 
     33   INTEGER , PARAMETER ::   jp_wspd = 7          ! index of 10 m neutral wind speed (m/s) 
     34 
     35   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wave ! structure of input fields (file informations, fields read) 
     36   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave ! CCC moved to sbc_oce 
     37   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wavepar ! structure of input fields (file informations, fields read) 
     38   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: ust_wavepar, vst_wavepar, swh_wavepar, mwp_wavepar 
     39   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: phioc_wavepar  !WAVE2NEMO 14.12.2016 moved rn_crban to sbc_oce 
     40   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: tauoc_wavepar ! CCC moved to sbc_oce 
     41   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: wspd_wavepar  ! CCC moved to sbc_oce 
     42   !REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdww_wavepar 
    3643 
    3744   !! * Substitutions 
    3845#  include "domzgr_substitute.h90" 
     46#  include "vectopt_loop_substitute.h90" 
    3947   !!---------------------------------------------------------------------- 
    4048   !! NEMO/OPA 4.0 , NEMO Consortium (2011)  
     
    5159      !! 
    5260      !! ** Method  : - Read namelist namsbc_wave 
    53       !!              - Read Cd_n10 fields in netcdf files  
    54       !!              - Read stokes drift 2d in netcdf files  
    55       !!              - Read wave number      in netcdf files  
    56       !!              - Compute 3d stokes drift using monochromatic 
    57       !! ** action  :    
    58       !!                
    59       !!--------------------------------------------------------------------- 
    60       USE oce,  ONLY : un,vn,hdivn,rotn 
    61       USE divcur 
    62       USE wrk_nemo 
    63 #if defined key_bdy 
    64       USE bdy_oce, ONLY : bdytmask 
    65 #endif 
    66       INTEGER, INTENT( in  ) ::  kt       ! ocean time step 
     61      !!              - Read Cd_n10 fields in netcdf files 
     62      !! ** action  : 
     63      !! 
     64      !!--------------------------------------------------------------------- 
     65      INTEGER, INTENT( in  ) ::  kt   ! ocean time step 
    6766      INTEGER                ::  ierror   ! return error code 
    68       INTEGER                ::  ifpr, jj,ji,jk  
    69       INTEGER                ::   ios     ! Local integer output status for namelist read 
    70       REAL(wp),DIMENSION(:,:,:),POINTER             ::  udummy,vdummy,hdivdummy,rotdummy 
    71       REAL                                          ::  z2dt,z1_2dt 
    72       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i     ! array of namelist informations on the fields to read 
    73       CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files 
    74       TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_wn   ! informations about the fields to be read 
    75       !!--------------------------------------------------------------------- 
    76       NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_wn 
     67      CHARACTER(len=100)     ::  cn_dir_cdg                      ! Rootdirectory for location of drag coefficient files 
     68      TYPE(FLD_N)            ::  sn_cdg                          ! informationsabout the fields to be read 
     69      !!--------------------------------------------------------------------- 
     70      NAMELIST/namsbc_wave/  sn_cdg, cn_dir_cdg 
    7771      !!--------------------------------------------------------------------- 
    7872 
     
    8377      IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
    8478         !                                      ! -------------------- ! 
    85          REWIND( numnam_ref )              ! Namelist namsbc_wave in reference namelist : File for drag coeff. from wave model 
    86          READ  ( numnam_ref, namsbc_wave, IOSTAT = ios, ERR = 901) 
    87 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in reference namelist', lwp ) 
    88  
    89          REWIND( numnam_cfg )              ! Namelist namsbc_wave in configuration namelist : File for drag coeff. from wave model 
    90          READ  ( numnam_cfg, namsbc_wave, IOSTAT = ios, ERR = 902 ) 
    91 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave in configuration namelist', lwp ) 
    92          IF(lwm) WRITE ( numond, namsbc_wave ) 
    93          ! 
    94  
    95          IF ( ln_cdgw ) THEN 
    96             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    97             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    98             ! 
    99                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   ) 
    100             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) ) 
    101             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    102             ALLOCATE( cdn_wave(jpi,jpj) ) 
    103             cdn_wave(:,:) = 0.0 
    104         ENDIF 
    105          IF ( ln_sdw ) THEN 
    106             slf_i(jp_usd) = sn_usd ; slf_i(jp_vsd) = sn_vsd; slf_i(jp_wn) = sn_wn 
    107             ALLOCATE( sf_sd(3), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg 
    108             IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocate sf_wave structure' ) 
    109             ! 
    110             DO ifpr= 1, jpfld 
    111                ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) ) 
    112                IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) ) 
    113             END DO 
    114             CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' ) 
    115             ALLOCATE( usd2d(jpi,jpj),vsd2d(jpi,jpj),uwavenum(jpi,jpj),vwavenum(jpi,jpj) ) 
    116             ALLOCATE( usd3d(jpi,jpj,jpk),vsd3d(jpi,jpj,jpk),wsd3d(jpi,jpj,jpk) ) 
    117             usd2d(:,:) = 0.0 ;  vsd2d(:,:) = 0.0 ; uwavenum(:,:) = 0.0 ; vwavenum(:,:) = 0.0 
    118             usd3d(:,:,:) = 0.0 ;vsd3d(:,:,:) = 0.0 ; wsd3d(:,:,:) = 0.0 
     79         !                                            !* set file information 
     80         !                                            (default values) 
     81         ! ... default values (NB: frequency positive => hours, negative => 
     82         ! months) 
     83         !              !   file   ! frequency !  variable  ! time intep !  clim 
     84         !              ! 'yearly' or ! weights  ! rotation ! 
     85         !              !   name   !  (hours)  !   name     !   (T/F)    ! 
     86         !              (T/F)  !  'monthly'  ! filename ! pairs    ! 
     87         sn_cdg = FLD_N('cdg_wave'  ,    1     ,'drag_coeff',  .true.    ,.false. ,   'daily'   , ''       , '',''       ) 
     88         cn_dir_cdg = './'          ! directory in which the Patm data are 
     89 
     90 
     91         REWIND( numnam_ref )                             !* read in namlist namsbc_wave 
     92         READ  ( numnam_ref, namsbc_wave ) 
     93    REWIND( numnam_ref ) 
     94         REWIND( numnam_cfg )                             !* read in namlist namsbc_wave 
     95         READ  ( numnam_cfg, namsbc_wave ) 
     96    REWIND( numnam_cfg ) 
     97         ! 
     98 
     99         ALLOCATE( sf_wave(1), STAT=ierror )           !* allocate and fillsf_wave with sn_cdg 
     100         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave: unable to allocatesf_wave structure' ) 
     101         ! 
     102         CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave', 'Wavemodule ', 'namsbc_wave' ) 
     103                                ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1)   ) 
     104         IF( sn_cdg%ln_tint )   ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) ) 
     105         !ALLOCATE( cdn_wave(jpi,jpj) ) 
     106         ! Allocation done by sbc_oce 
     107         cdn_wave(:,:) = 0.0 
     108      ENDIF 
     109         ! 
     110         ! 
     111      CALL fld_read( kt, nn_fsbc, sf_wave )      !* read drag coefficient fromexternal forcing 
     112      cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1) 
     113   END SUBROUTINE sbc_wave 
     114 
     115 
     116   SUBROUTINE sbc_wavepar( kt ) 
     117      !!--------------------------------------------------------------------- 
     118      !!                     ***  ROUTINE sbc_wavepar  *** 
     119      !! 
     120      !! ** Purpose :   Provide at each time step wave model parameters, including the 
     121      !!                Stokes drift (east and north components), significant wave height and the 
     122      !!                mean wave period as well as the normalized stress and energy flux into the  
     123      !!                ocean for TKE. 
     124      !! 
     125      !! 
     126      !! ** Method  : - Read namelist namsbc_wavepar 
     127      !!              - Read fields in NetCDF files  
     128      !! ** action  :    
     129      !!                
     130      !!--------------------------------------------------------------------- 
     131      INTEGER, INTENT( in  ) ::  kt  ! ocean time step 
     132      !! 
     133      INTEGER  ::  ierror ! return error code 
     134      INTEGER  ::   ji,iii,jjj      ! dummy loop index, Victor added iii, jjj 
     135      INTEGER  ::   jfld    ! dummy loop arguments 
     136      !! 
     137      CHARACTER(len=100)     ::  cn_dir_wavepar ! Root directory for location of ECWAM wave parameter fields 
     138      TYPE(FLD_N), DIMENSION(jpfld_wavepar) :: slf_i ! array of namelist informations on the fields to read 
     139      TYPE(FLD_N)            ::  sn_ust, sn_vst, sn_swh, sn_mwp ! information about the fields to be read 
     140      TYPE(FLD_N)            ::  sn_phioc, sn_tauoc 
     141      TYPE(FLD_N)            ::  sn_wspd 
     142      !TYPE(FLD_N)            ::  sn_cdww 
     143      !!--------------------------------------------------------------------- 
     144      !NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_phioc, sn_tauoc, cn_dir_wavepar 
     145      NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_wspd, sn_phioc, sn_tauoc, cn_dir_wavepar 
     146      !!--------------------------------------------------------------------- 
     147 
     148      !!---------------------------------------------------------------------- 
     149      ! 
     150      ! 
     151      !                                         ! -------------------- ! 
     152      IF( kt == nit000 ) THEN                   ! First call kt=nit000 ! 
     153         !                                      ! -------------------- ! 
     154         ! set file information (default values) 
     155         ! ... default values (NB: frequency positive => hours, negative => months) 
     156         !            !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation ! 
     157         !            !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    ! 
     158         sn_ust = FLD_N( 'ust'   ,     6     ,  'ust'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''       ) 
     159         sn_vst = FLD_N( 'vst'   ,     6     ,  'vst'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''      ) 
     160         sn_swh = FLD_N( 'swh'   ,     6     ,  'swh'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''      ) 
     161         sn_mwp = FLD_N( 'mwp'   ,     6     ,  'mwp'     ,  .true.    , .false. ,   'monthly' , ''       , '' ,''      ) 
     162         !sn_cdww = FLD_N( 'cdww'  ,    6     ,  'cdww'    ,  .true.    , .false. ,   'monthly' , ''       , ''       ) 
     163         sn_wspd = FLD_N( 'wspd'   ,   6     ,  'wspd'    ,  .true.    , .false. ,   'monthly' , ''       , '',''       ) 
     164         sn_phioc = FLD_N( 'phioc'   , 6     ,  'phioc'   ,  .true.    , .false. ,   'monthly' , ''       , '',''       ) 
     165         sn_tauoc = FLD_N( 'tauoc'   , 6     ,  'tauoc'   ,  .true.    , .false. ,   'monthly' , ''       , '',''       ) 
     166         cn_dir_wavepar = './'          ! directory in which the wave data are found 
     167 
     168         REWIND( numnam_ref )                             !* read in namlist namsbc_wave 
     169         READ  ( numnam_ref, namsbc_wavepar )  
     170         REWIND ( numnam_ref ) 
     171 
     172         REWIND( numnam_cfg )                             !* read in namlist namsbc_wave 
     173         READ  ( numnam_cfg, namsbc_wavepar )  
     174         REWIND ( numnam_cfg ) 
     175         ! 
     176 
     177 
     178         slf_i(jp_ust) = sn_ust  
     179         slf_i(jp_vst) = sn_vst 
     180         slf_i(jp_swh) = sn_swh 
     181         slf_i(jp_mwp) = sn_mwp 
     182         !slf_i(jp_cdww) = sn_cdww 
     183         slf_i(jp_wspd) = sn_wspd 
     184         slf_i(jp_phioc) = sn_phioc 
     185         slf_i(jp_tauoc) = sn_tauoc 
     186 
     187         ALLOCATE( sf_wavepar(jpfld_wavepar), STAT=ierror )          !* set sf_wavepar structure 
     188         IF ( ierror > 0 ) THEN 
     189            CALL ctl_stop( 'STOP', 'sbc_wavepar: unable to allocate sf_wavepar structure' ) ; RETURN 
    119190         ENDIF 
     191         ! 
     192         jfld = jpfld_wavepar 
     193         DO ji = 1, jpfld_wavepar 
     194            ALLOCATE( sf_wavepar(ji)%fnow(jpi,jpj,1)   ) 
     195            IF ( slf_i(ji)%ln_tint ) ALLOCATE( sf_wavepar(ji)%fdta(jpi,jpj,1,2) ) 
     196         ENDDO 
     197         ALLOCATE( ust_wavepar(jpi,jpj) ) 
     198         ALLOCATE( vst_wavepar(jpi,jpj) ) 
     199         ALLOCATE( swh_wavepar(jpi,jpj) ) 
     200         ALLOCATE( mwp_wavepar(jpi,jpj) ) 
     201         !ALLOCATE( cdww_wavepar(jpi,jpj) ) 
     202         !ALLOCATE( wspd_wavepar(jpi,jpj) ) 
     203         ALLOCATE( phioc_wavepar(jpi,jpj) ) 
     204         ALLOCATE( tauoc_wavepar(jpi,jpj) ) 
     205         ! 
     206         CALL fld_fill( sf_wavepar, slf_i, cn_dir_wavepar, 'sbc_wavepar', 'Wave module ', 'namsbc_wavepar' ) 
     207         ! 
    120208      ENDIF 
    121209         ! 
    122          ! 
    123       IF ( ln_cdgw ) THEN 
    124          CALL fld_read( kt, nn_fsbc, sf_cd )      !* read drag coefficient from external forcing 
    125          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1) 
    126       ENDIF 
    127       IF ( ln_sdw )  THEN 
    128           CALL fld_read( kt, nn_fsbc, sf_sd )      !* read drag coefficient from external forcing 
    129  
    130          ! Interpolate wavenumber, stokes drift into the grid_V and grid_V 
    131          !------------------------------------------------- 
    132  
    133          DO jj = 1, jpjm1 
    134             DO ji = 1, jpim1 
    135                uwavenum(ji,jj)=0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    136                &                                + sf_sd(3)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    137  
    138                vwavenum(ji,jj)=0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(3)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    139                &                                + sf_sd(3)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    140  
    141                usd2d(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( sf_sd(1)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    142                &                                + sf_sd(1)%fnow(ji+1,jj,1) * tmask(ji+1,jj,1) ) 
    143  
    144                vsd2d(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( sf_sd(2)%fnow(ji,jj,1) * tmask(ji,jj,1) & 
    145                &                                + sf_sd(2)%fnow(ji,jj+1,1) * tmask(ji,jj+1,1) ) 
    146             END DO 
    147          END DO 
    148  
    149           !Computation of the 3d Stokes Drift 
    150           DO jk = 1, jpk 
    151              DO jj = 1, jpj-1 
    152                 DO ji = 1, jpi-1 
    153                    usd3d(ji,jj,jk) = usd2d(ji,jj)*exp(2.0*uwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk)))) 
    154                    vsd3d(ji,jj,jk) = vsd2d(ji,jj)*exp(2.0*vwavenum(ji,jj)*(-MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk)))) 
    155                 END DO 
    156              END DO 
    157              usd3d(jpi,:,jk) = usd2d(jpi,:)*exp( 2.0*uwavenum(jpi,:)*(-gdept_0(jpi,:,jk)) ) 
    158              vsd3d(:,jpj,jk) = vsd2d(:,jpj)*exp( 2.0*vwavenum(:,jpj)*(-gdept_0(:,jpj,jk)) ) 
    159           END DO 
    160  
    161           CALL wrk_alloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    162            
    163           udummy(:,:,:)=un(:,:,:) 
    164           vdummy(:,:,:)=vn(:,:,:) 
    165           hdivdummy(:,:,:)=hdivn(:,:,:) 
    166           rotdummy(:,:,:)=rotn(:,:,:) 
    167           un(:,:,:)=usd3d(:,:,:) 
    168           vn(:,:,:)=vsd3d(:,:,:) 
    169           CALL div_cur(kt) 
    170       !                                           !------------------------------! 
    171       !                                           !     Now Vertical Velocity    ! 
    172       !                                           !------------------------------! 
    173           z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    174  
    175           z1_2dt = 1.e0 / z2dt 
    176           DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence 
    177              ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 
    178              wsd3d(:,:,jk) = wsd3d(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        & 
    179                 &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    & 
    180                 &                         * tmask(:,:,jk) * z1_2dt 
    181 #if defined key_bdy 
    182              wsd3d(:,:,jk) = wsd3d(:,:,jk) * bdytmask(:,:) 
    183 #endif 
    184           END DO 
    185           hdivn(:,:,:)=hdivdummy(:,:,:) 
    186           rotn(:,:,:)=rotdummy(:,:,:) 
    187           vn(:,:,:)=vdummy(:,:,:) 
    188           un(:,:,:)=udummy(:,:,:) 
    189           CALL wrk_dealloc( jpi,jpj,jpk,udummy,vdummy,hdivdummy,rotdummy) 
    190       ENDIF 
    191    END SUBROUTINE sbc_wave 
     210      CALL fld_read( kt, nn_fsbc, sf_wavepar )      !* read wave parameters from ECWAM NetCDF file 
     211      ust_wavepar(:,:) = sf_wavepar(jp_ust)%fnow(:,:,1) 
     212      vst_wavepar(:,:) = sf_wavepar(jp_vst)%fnow(:,:,1) 
     213      swh_wavepar(:,:) = sf_wavepar(jp_swh)%fnow(:,:,1) 
     214      mwp_wavepar(:,:) = sf_wavepar(jp_mwp)%fnow(:,:,1) 
     215      !cdww_wavepar(:,:) = sf_wavepar(jp_cdww)%fnow(:,:,1) 
     216      wspd_wavepar(:,:) = sf_wavepar(jp_wspd)%fnow(:,:,1) 
     217      phioc_wavepar(:,:) = sf_wavepar(jp_phioc)%fnow(:,:,1) 
     218      tauoc_wavepar(:,:) = sf_wavepar(jp_tauoc)%fnow(:,:,1) 
     219  
     220      
     221       rn_crban(:,:)=29*phioc_wavepar(:,:) ! Alfa is phioc*sqrt(rw/ra)sbc_wa 
     222       ! Limit cr_ban between 10 and 300 
     223       DO iii=1,jpi 
     224           DO jjj=1,jpj 
     225 
     226                IF (rn_crban(iii,jjj) < 10 ) THEN 
     227                rn_crban(iii,jjj)=10                 
     228                ELSEIF (rn_crban(iii,jjj) > 300) THEN                
     229                rn_crban(iii,jjj)=300 
     230                ENDIF 
     231         
     232           ENDDO 
     233       ENDDO 
     234  
     235     END SUBROUTINE sbc_wavepar 
    192236       
    193    !!====================================================================== 
     237 
    194238END MODULE sbcwave 
  • 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  
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6204 r7807  
    323323                               CALL dyn_vor( kstp )         ! vorticity term including Coriolis 
    324324                               CALL dyn_ldf( kstp )         ! lateral mixing 
     325        IF ( ln_stcor )        CALL dyn_stcor ( kstp )      ! WAVE2NEMO add Stokes-Coriolis forcing 
    325326        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified) 
    326327#if defined key_agrif 
  • branches/UKMO/r6232_HZG_WAVE/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r6204 r7807  
    110110   USE timing           ! Timing 
    111111 
     112   USE dynstcor         ! WAVE2NEMO Stokes-Coriolis forcing 
    112113#if defined key_agrif 
    113114   USE agrif_opa_sponge ! Momemtum and tracers sponges 
Note: See TracChangeset for help on using the changeset viewer.