Changeset 4351


Ignore:
Timestamp:
2014-01-16T15:27:29+01:00 (7 years ago)
Author:
acc
Message:

Branch 2013/dev_ECMWF_waves ticket #1216. First set of code imports from ECMWF. Untested

Location:
branches/2013/dev_ECMWF_waves/NEMOGCM
Files:
2 added
9 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_ECMWF_waves/NEMOGCM/ARCH/arch-ALTIX_NAUTILUS_MPT.fcm

    r4306 r4351  
    4747%LD                  ifort 
    4848%FPPFLAGS            -P -C -traditional 
    49 %LDFLAGS             -lmpi -lstdc++ -lcurl 
     49%LDFLAGS             -lmpi -lstdc++ -lcurl -L/fibre/acc/PIN2CPU -lsgipin 
    5050%AR                  ar  
    5151%ARFLAGS             -r 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg

    r4257 r4351  
    184184!----------------------------------------------------------------------- 
    185185/ 
     186!----------------------------------------------------------------------- 
     187&namsbc_wave  ! External fields from wave model 
     188!----------------------------------------------------------------------- 
     189/ 
     190!----------------------------------------------------------------------- 
     191&namsbc_wave_ecmwf   ! ECWMF external wave model drag coefficient 
     192!----------------------------------------------------------------------- 
     193/ 
     194!----------------------------------------------------------------------- 
     195&namsbc_wavepar ! namsbc_wavepar  parameters from wave model 
     196!----------------------------------------------------------------------- 
     197/ 
     198!----------------------------------------------------------------------- 
     199&namsbc_waveparlim 
     200!----------------------------------------------------------------------- 
     201/ 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4347 r4351  
    225225   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave) 
    226226   ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave) 
     227   ln_stcor    = .false.   !   Stokes drift read from wave model 
     228   ln_wavetke  = .false.   !   Wave parameters from wave model for the TKE BC 
     229   ln_tauoc    = .false.   !   Wave-modified stress from wave model 
    227230   cn_iceflx = 'linear'    !  redistribution of solar input into ice categories during coupling ice/atm. 
    228231/ 
     
    11461149/ 
    11471150!----------------------------------------------------------------------- 
     1151&namsbc_wave_ecmwf   ! External fields from wave model 
     1152!----------------------------------------------------------------------- 
     1153!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     1154!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
     1155   sn_cdg      =  'cdww' ,        6          , 'cdww'       ,     .true.   , .false. , 'monthly'  , ''      , ''       , '' ! 
     1156   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files 
     1157/ 
     1158!----------------------------------------------------------------------- 
    11481159&namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed) 
    11491160!----------------------------------------------------------------------- 
     
    11601171   rn_htrmax         =  200.0   ! max. depth of transition range 
    11611172/ 
     1173!----------------------------------------------------------------------- 
     1174&namsbc_wavepar ! namsbc_wavepar  parameters from wave model (CCC not identical to namsbc_wave) 
     1175!----------------------------------------------------------------------- 
     1176!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! 
     1177!                                !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! 
     1178   sn_ust      = 'ust'            ,         6         , 'ust'     ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1179   sn_vst      = 'vst'            ,         6         , 'vst'     ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1180   sn_swh      = 'swh'            ,         6         , 'swh'     ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1181   sn_mwp      = 'mwp'            ,         6         , 'mwp'     ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1182   sn_wspd     = 'wspd'           ,         6         , 'wspd'    ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1183   sn_phioc    = 'physphioc'      ,         6         , 'physphioc',  .true.    , .false. , 'monthly'  , ''       , '' 
     1184   sn_tauoc    = 'tauoc'          ,         6         , 'tauoc'   ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1185   cn_dir_wavepar  = './'      !  root directory for the location of the ECWAM files 
     1186/ 
     1187!----------------------------------------------------------------------- 
     1188&namsbc_waveparlim 
     1189!----------------------------------------------------------------------- 
     1190!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! 
     1191!                                !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! 
     1192   sn_strn   = 'istr'            ,         6         , 'icestrain'     ,   .true.    , .false. , 'monthly'  , ''       , '' 
     1193   cn_dir_waveparlim  = './'      !  root directory for the location of the ECWAM files 
     1194/ 
     1195!----------------------------------------------------------------------- 
     1196&namsbc_wave   ! External fields from wave model 
     1197!----------------------------------------------------------------------- 
     1198sn_cdg     = 'cdww' ,        6          , 'cdww' , .true.   , .false. , 'monthly'  ,''         , ''  ! 
     1199   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files 
     1200/ 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r4306 r4351  
    5353   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model 
    5454   ! 
     55   REAL(wp), PUBLIC ::   rn_rrelwind = 1.0       !: relative wind ratio (1.0 means surface currents are subtracted, 0.0 means no current) 
     56   LOGICAL , PUBLIC ::   ln_wavetke  = .FALSE.   !: true if wave parameters are read 
     57   LOGICAL , PUBLIC ::   ln_stcor    = .FALSE.   !: true if Stokes-Coriolis forcing is included  
     58   LOGICAL , PUBLIC ::   ln_tauoc    = .FALSE.   !: true if wave-modified water-side stress is used 
     59 
    5560   LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs 
    5661   ! 
     
    6974   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2] 
    7075   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2]  
     76   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_wave          !: wave model drag coefficient [N/m2]  
     77   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wspd_wavepar      !: wave model 10-m neutral wind [m/s] CCC 
     78   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tauoc_wavepar     !: normalized stress into the ocean from wave model  CCC 
    7179   !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads 
    7280   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s] 
     
    131139      ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
    132140         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
     141         ! Initialize these since that is not done in the code. 
     142      rnf        (:,:) = 0.0_wp 
     143      sbc_tsc  (:,:,:) = 0.0_wp 
     144      qsr_hc   (:,:,:) = 0.0_wp 
     145      rnf_b      (:,:) = 0.0_wp 
     146      sbc_tsc_b(:,:,:) = 0.0_wp 
     147      qsr_hc_b (:,:,:) = 0.0_wp 
     148 
    133149         ! 
    134150      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
     
    138154         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       & 
    139155         &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) ) 
     156      ALLOCATE( cdn_wave(jpi,jpj) ) 
     157      ALLOCATE( wspd_wavepar(jpi,jpj) ) 
    140158         ! 
    141159#if defined key_vvl 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r4333 r4351  
    3838   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3939   USE prtctl          ! Print control 
    40    USE sbcwave,ONLY :  cdn_wave !wave module  
     40   !USE sbcwave,ONLY :  cdn_wave !wave module  ! moved to sbc_oce  
    4141#if defined key_lim3 || defined key_cice 
    4242   USE sbc_ice         ! Surface boundary condition: ice fields 
     
    6464    
    6565   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read) 
     66   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Cd_n10 ! 10m neutral drag coefficient for output 
    6667          
    6768   !                                             !!! CORE bulk parameters 
     
    8384   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    8485   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     86   LOGICAL  ::   ln_cdec   = .FALSE.   ! logical flag for using ec neutral wind drag coef  
     87   ! References : P. Janssen, 2008, ECMWF Workshop on Ocean-Atmosphere Interactions, 10-12 November p47-60 
    8588 
    8689   !! * Substitutions 
     
    138141      NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  & 
    139142         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
     143         &                  sn_ccov, ln_cdec, ln_cldcov, rn_rrelwind 
    140144         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    141145         &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu 
     
    163167            sn_qsr%ln_tint = .false. 
    164168         ENDIF 
     169         IF( ln_cdec .AND. lwp ) WRITE(numout,*)'Using Hans Hersbach formula for drag.' 
    165170         !                                         ! store namelist information in an array 
    166171         slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
     
    182187         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' ) 
    183188         ! 
     189         ! Drag coefficent so we can write to disk 
     190         ALLOCATE(Cd_n10(jpi,jpj)) 
     191 
    184192         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90) 
    185193         ! 
     
    239247      INTEGER  ::   ji, jj               ! dummy loop indices 
    240248      REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable 
     249      REAL(wp) ::   ztheta               ! local variable, wind direction 
    241250      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    242251      REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst 
     
    283292!CDIR COLLAPSE 
    284293#endif 
    285       DO jj = 2, jpjm1 
    286          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    287             zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    288             zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     294      ! Wind vector at T points 
     295      IF ( ln_cdgw ) THEN 
     296         ! Use neutral 10-m wind speed from wave model if available 
     297         DO jj = 2, jpjm1 
     298            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     299               ! Local wind direction from non-neutral 10-m wind vector (relative to grid and no current) 
     300               ! since we do not have the directional information from the wave model 
     301               ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1), sf(jp_wndj)%fnow(ji,jj,1)) 
     302               ! Wind vector magnitude from 10-m neutral wind speed from wave model 
     303               zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta) 
     304               zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta) 
     305               ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0 
     306               zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_rrelwind*(pu(ji-1,jj  ) + pu(ji,jj))) 
     307               zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_rrelwind*(pv(ji  ,jj-1) + pv(ji,jj))) 
     308            END DO 
    289309         END DO 
    290       END DO 
     310      ELSE 
     311         DO jj = 2, jpjm1 
     312            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     313               ! Correct for surface current, 0.0 <= rn_rrelwind <= 1.0 
     314               zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - 0.5*rn_rrelwind*( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
     315               zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - 0.5*rn_rrelwind*( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
     316            END DO 
     317         END DO 
     318      ENDIF 
    291319      CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. ) 
    292320      CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. ) 
     
    353381 
    354382      ! ... tau module, i and j component 
    355       DO jj = 1, jpj 
    356          DO ji = 1, jpi 
    357             zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 
    358             taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    359             zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    360             zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
    361          END DO 
    362       END DO 
     383      ! Use cdn_wave instead of Cd if we have wave parameters 
     384      IF ( ln_cdgw ) THEN 
     385         DO jj = 1, jpj 
     386            DO ji = 1, jpi 
     387               !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj) 
     388               zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj) 
     389               !taum(ji,jj) = tauoc_wavepar(ji,jj) 
     390               !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj) 
     391               taum  (ji,jj) = zztmp * wndm  (ji,jj) 
     392               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     393               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     394            ENDDO 
     395         ENDDO 
     396      ELSE 
     397         DO jj = 1, jpj 
     398            DO ji = 1, jpi 
     399               zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj) 
     400               taum  (ji,jj) = zztmp * wndm  (ji,jj) 
     401               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     402               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     403            ENDDO 
     404         ENDDO 
     405      ENDIF 
    363406 
    364407      ! ... add the HF tau contribution to the wind stress module? 
     
    379422      CALL lbc_lnk( utau(:,:), 'U', -1. ) 
    380423      CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
     424 
     425      IF (ln_tauoc) THEN  
     426         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:) 
     427         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:) 
     428         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:) 
     429      ENDIF 
    381430 
    382431      !  Turbulent fluxes over ocean 
     
    755804      REAL(wp), DIMENSION(:,:), POINTER  ::   dT            ! air/sea temperature differeence      [K] 
    756805      REAL(wp), DIMENSION(:,:), POINTER  ::   dq            ! air/sea humidity difference          [K] 
    757       REAL(wp), DIMENSION(:,:), POINTER  ::   Cd_n10        ! 10m neutral drag coefficient 
    758806      REAL(wp), DIMENSION(:,:), POINTER  ::   Ce_n10        ! 10m neutral latent coefficient 
    759807      REAL(wp), DIMENSION(:,:), POINTER  ::   Ch_n10        ! 10m neutral sensible coefficient 
     
    775823      ! 
    776824      CALL wrk_alloc( jpi,jpj, stab )   ! integer 
    777       CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
     825      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    778826      CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 
    779827 
    780828      !! * Start 
    781829      !! Air/sea differences 
    782       dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s 
     830      dU10 = MAX(0.5, dU)     ! we don't want to fall under 0.5 m/s 
    783831      dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu 
    784832      dq = q_a - q_sat 
     
    788836      !! 
    789837      !! Neutral Drag Coefficient 
    790       stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0  
    791       IF  ( ln_cdgw ) THEN 
    792         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 
    793         Cd_n10(:,:) =   cdn_wave 
     838      stab   = 0.5 + SIGN(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0  
     839 
     840      ! Use drag coefficient from a wave model ... 
     841      IF (ln_cdgw) THEN 
     842         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 
     843         Cd_n10(:,:) =  cdn_wave(:,:) 
    794844      ELSE 
    795         Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a) 
     845      ! ... or choose a drag law 
     846         IF (ln_cdec) THEN 
     847            Cd_n10  = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21)    ! PJ (6) 
     848         ELSE 
     849            Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )              ! L & Y eq. (6a) 
     850         ENDIF 
    796851      ENDIF 
    797852      sqrt_Cd_n10 = sqrt(Cd_n10) 
     
    824879 
    825880           !! Updating the neutral 10m transfer coefficients : 
    826            Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a) 
     881           IF (ln_cdec) THEN 
     882             Cd_n10  = 1E-3 * ( 1.03 + 0.04 * U_n10**1.48 ) / (U_n10**0.21)    ! PJ (6) 
     883           ELSE 
     884             Cd_n10  = 1E-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 )              ! L & Y eq. (6a) 
     885           ENDIF 
    827886           sqrt_Cd_n10 = sqrt(Cd_n10) 
    828887           Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b) 
     
    847906      !! 
    848907      CALL wrk_dealloc( jpi,jpj, stab )   ! integer 
    849       CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
     908      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    850909      CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m ) 
    851910      ! 
     
    891950      REAL(wp), DIMENSION(:,:), POINTER ::   dT            ! air/sea temperature differeence   [K] 
    892951      REAL(wp), DIMENSION(:,:), POINTER ::   dq            ! air/sea humidity difference       [K] 
    893       REAL(wp), DIMENSION(:,:), POINTER ::   Cd_n10        ! 10m neutral drag coefficient 
    894952      REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient 
    895953      REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient 
     
    911969      IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_2Z') 
    912970      ! 
    913       CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
     971      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    914972      CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 
    915973      CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) 
     
    923981      !! Neutral Drag Coefficient : 
    924982      stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE 
    925       IF( ln_cdgw ) THEN 
    926         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 
    927         Cd_n10(:,:) =   cdn_wave 
     983 
     984      ! Use drag coefficient from a wave model ... 
     985      IF (ln_cdgw) THEN 
     986         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1) 
     987         Cd_n10(:,:) =  cdn_wave 
    928988      ELSE 
    929         Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 )  
     989      ! ... or choose a drag law 
     990         IF (ln_cdec) THEN 
     991            Cd_n10  = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21)    ! PJ (6) 
     992         ELSE 
     993            Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )              ! L & Y eq. (6a) 
     994         ENDIF 
    930995      ENDIF 
    931996      sqrt_Cd_n10 = sqrt(Cd_n10) 
     
    9961061      END DO 
    9971062      !! 
    998       CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
     1063      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L ) 
    9991064      CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 ) 
    10001065      CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m ) 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4333 r4351  
    5050   USE timing           ! Timing 
    5151   USE sbcwave          ! Wave module 
     52   USE sbcwave_ecmwf    ! ECMWF wave module 
     53   USE sbcwavelim       ! Wave for LIM module 
    5254 
    5355   IMPLICIT NONE 
     
    8284      INTEGER ::   icpt   ! local integer 
    8385      !! 
    84       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   & 
    85          &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   & 
    86          &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, cn_iceflx 
     86      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx ,  ln_blk_clio, ln_blk_core, ln_cpl,   & 
     87         &             ln_blk_mfs, ln_apr_dyn, nn_ice ,  nn_ice_embd, ln_dm2dc   , ln_rnf,   & 
     88         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave    , ln_sdw     , nn_lsm,   & 
     89         &             cn_iceflx , ln_stcor  , ln_tauoc, ln_wavetke  
    8790      INTEGER  ::   ios 
    8891      !!---------------------------------------------------------------------- 
     
    135138         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea 
    136139         WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm 
     140         WRITE(numout,*) '              Drag coefficient from wave model           ln_cdgw     = ', ln_cdgw   
     141         WRITE(numout,*) '              Stokes-Coriolis forcing from wave model    ln_stcor    = ', ln_stcor   
     142         WRITE(numout,*) '              TKE wave breaking BC from wave model       ln_wavetke  = ', ln_wavetke 
     143         WRITE(numout,*) '              Wave-modified stress from wave model       ln_tauoc    = ', ln_tauoc 
     144         WRITE(numout,*) '              Wave-modified stress min. sea. frac     rn_taouc_cimin = ', rn_taouc_cimin 
    137145      ENDIF 
    138146 
     
    216224      ELSE 
    217225      IF ( ln_cdgw .OR. ln_sdw  )                                         &  
    218          &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     & 
     226         &   CALL ctl_warn('Not Activated Wave Module (ln_wave=F) but     & 
    219227         & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ') 
    220228      ENDIF  
     
    300308      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m) 
    301309      !                                                  ! averaged over nf_sbc time-step 
     310 
     311      IF (ln_cdgw) CALL sbc_wave_ecmwf( kt ) 
     312      ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for TKE are needed 
     313      ! Also read it if drag coefficient to ensure that we use the neutral 10m from WAM. 
     314      IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc .OR. ln_cdgw) THEN 
     315         CALL sbc_wavepar( kt ) 
     316      ENDIF 
    302317 
    303318      IF (ln_wave) CALL sbc_wave( kt ) 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r4147 r4351  
    2626   !!                 !                                + cleaning of the parameters + bugs correction 
    2727   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     28   !!                 !  2012-12  (oyvind.breivik@ecwmf.int) Changed to wave breaking source term 
    2829   !!---------------------------------------------------------------------- 
    2930#if defined key_zdftke   ||   defined key_esopa 
     
    4243   USE domvvl         ! domain: variable volume layer 
    4344   USE sbc_oce        ! surface boundary condition: ocean 
     45   USE sbcwave        ! wave parameters 
     46   USE sbcwave_ecmwf  ! ECMWF wave parameters 
    4447   USE zdf_oce        ! vertical physics: ocean variables 
    4548   USE zdfmxl         ! vertical physics: mixed layer 
     
    7881   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not 
    7982   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells 
     83   REAL(wp) ::   rn_alpavg = 100.0_wp      ! Average alpha (normalized, water-side energy flux, only used in 
     84                                           ! case we use the source term formulation with no wave model information 
     85                                           ! NB! The option of running with source term and no wave model has 
     86                                           ! not been properly tested yet but is included for consistency (Oyvind Breivik, 2013-07-16) 
     87   REAL(wp) ::   rn_deptmaxtkew = 50.0_wp   ! maximum depth [m] to be affected by TKE from breaking waves  
    8088 
    8189   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 
     
    8593 
    8694   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2] 
     95   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   e0             !: top level turbulent kinetic energy   [m2/s2] from waves 
     96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dfac           !: wave depth scaling factor    
    8797   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau) 
    8898   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation 
     
    109119      !!                ***  FUNCTION zdf_tke_alloc  *** 
    110120      !!---------------------------------------------------------------------- 
     121      INTEGER :: iexalloc 
    111122      ALLOCATE(                                                                    & 
    112123#if defined key_c1d 
     
    118129         &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      ) 
    119130         ! 
     131      IF ( ln_wavetke ) THEN 
     132         ALLOCATE( e0(jpi,jpj) , dfac(jpi,jpj), STAT= iexalloc ) 
     133         ! Setting arrays to zero initially. 
     134         e0(:,:) = 0.0_wp 
     135         dfac(:,:) = 0.0_wp 
     136         ! 
     137         zdf_tke_alloc = zdf_tke_alloc + iexalloc 
     138      ENDIF 
    120139      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc ) 
    121140      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays') 
     
    219238      REAL(wp) ::   zus   , zwlc  , zind            !    -         - 
    220239      REAL(wp) ::   zzd_up, zzd_lw                  !    -         - 
     240      REAL(wp) ::   lambda !    -         - 
     241      REAL(wp) ::   zphio_flux              !    -         - 
     242      REAL(wp) ::   sbr, z0, fb                     !    -         - 
     243 
    221244!!bfr      REAL(wp) ::   zebot                           !    -         - 
    222245      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc 
     
    239262      !                     !  Surface boundary condition on tke 
    240263      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    241       DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242          DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    244          END DO 
    245       END DO 
    246        
     264      ! Wave parameters read? 
     265      IF ( ln_wavetke ) THEN 
     266         ! Breaking-wave TKE injection as flux 
     267         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
     268            DO ji = fs_2, fs_jpim1   ! vector opt. 
     269               ! Roughness length proportional to wave height, capped at Hs = 0.1 m 
     270               z0 = 0.5_wp*MAX(swh_wavepar(ji,jj), 0.1_wp) 
     271               lambda = 2.38_wp/z0 
     272               ! The energy flux from the wave model in physical units 
     273               zphio_flux = phioc_wavepar(ji,jj) 
     274               ! TKE at surface (Mellor and Blumberg, 2004) 
     275               e0(ji,jj) = 0.5_wp*( 15.8_wp*zphio_flux/rau0 )**0.67_wp 
     276               ! Reducing the surface TKE by scaling with T-depth of first vertical level 
     277               dfac(ji,jj) = MAX( (2.0_wp*lambda*fsdept(ji,jj,1) )/3.0_wp, 0.00001_wp ) 
     278               ! TKE surface boundary condition weighted by the thickness of the first level 
     279               en(ji,jj,1) = MAX( rn_emin0, e0(ji,jj) * ( 1.0_wp - EXP(-dfac(ji,jj)) ) / dfac(ji,jj) ) * tmask(ji,jj,1) 
     280            ENDDO 
     281         ENDDO 
     282         ! If wave parameters are not read, use rn_ebb, corresponding to alpha = 100.0 
     283      ELSE 
     284         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
     285            DO ji = fs_2, fs_jpim1   ! vector opt. 
     286               ! Estimate average energy flux from std value found in rn_ebb (corresponds to alpha=100.0) 
     287               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     288            ENDDO 
     289         ENDDO 
     290      ENDIF 
     291  
    247292!!bfr   - start commented area 
    248293      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     
    443488      ENDIF 
    444489      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
     490      IF ( ln_wavetke ) THEN 
     491         CALL lbc_lnk( e0,     'W', 1. ) ! Lateral boundary conditions (sign unchanged) 
     492         CALL lbc_lnk( dfac,   'W', 1. ) ! Lateral boundary conditions (sign unchanged) 
     493      ENDIF 
    445494      ! 
    446495      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4338 r4351  
    284284                               CALL dyn_ldf( kstp )         ! lateral mixing 
    285285        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified) 
     286      IF( ln_stcor )           CALL dyn_stcor ( kstp )      ! Stokes-Coriolis forcing (ln_stcor set in SBC/sbc_oce.F90) 
    286287#if defined key_agrif 
    287288        IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge 
  • branches/2013/dev_ECMWF_waves/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r4328 r4351  
    5151   USE dynspg           ! surface pressure gradient        (dyn_spg routine) 
    5252   USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine) 
     53   USE dynstcor         ! Stokes-Coriolis forcing          (dyn_stcor routine) 
    5354 
    5455   USE dynnxt           ! time-stepping                    (dyn_nxt routine) 
Note: See TracChangeset for help on using the changeset viewer.