Changeset 10478

Timestamp:

2019-01-09T12:11:55+01:00 (5 years ago)

Author:

jcastill

Message:

Merged branch UKMO/r6232_HZG_WAVE-coupling@9868

Location:

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM

Files:

• CONFIG/SHARED/namelist_ref (modified) (9 diffs)
• NEMO/OPA_SRC/DIA/diawri.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/DYN/dynvor.F90 (modified) (1 diff)
• NEMO/OPA_SRC/SBC/sbc_oce.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SBC/sbcapr.F90 (modified) (5 diffs)
• NEMO/OPA_SRC/SBC/sbcblk_core.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (25 diffs)
• NEMO/OPA_SRC/SBC/sbcflx.F90 (modified) (5 diffs)
• NEMO/OPA_SRC/SBC/sbcmod.F90 (modified) (8 diffs)
• NEMO/OPA_SRC/SBC/sbcssm.F90 (modified) (9 diffs)
• NEMO/OPA_SRC/SBC/sbcwave.F90 (modified) (14 diffs)
• NEMO/OPA_SRC/ZDF/zdf_oce.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/ZDF/zdfgls.F90 (modified) (18 diffs)
• NEMO/OPA_SRC/ZDF/zdfini.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/step.F90 (modified) (1 diff)

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/CONFIG/SHARED/namelist_ref

@@ -240 +240 @@
    ln_blk_core = .true.    !  CORE bulk formulation                     (T => fill namsbc_core)
    ln_blk_mfs  = .false.   !  MFS bulk formulation                      (T => fill namsbc_mfs )
-   ln_cpl      = .false.   !  atmosphere coupled   formulation          ( requires key_oasis3 )
-   ln_mixcpl   = .false.   !  forced-coupled mixed formulation          ( requires key_oasis3 )
+   ln_cpl      = .false.   !  coupled   formulation                       ( requires key_oasis3 ) 
+   ln_mixcpl   = .false.   !  forced-coupled mixed atmosphere formulation ( requires key_oasis3 ) 
+   ln_wavcpl   = .false.   !  forced-coupled mixed wave       formulation ( requires key_oasis3 )
    nn_components = 0       !  configuration of the opa-sas OASIS coupling
                            !  =0 no opa-sas OASIS coupling: default single executable configuration
@@ -265 +266 @@
                            !     =2 annual global mean of e-p-r set to zero
    ln_wave = .false.       !  Activate coupling with wave (T => fill namsbc_wave)   
-   ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)   
-   ln_sdw  = .false.       !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)    
-   ln_tauoc= .false.       !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)   
-   ln_stcor= .false.       !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave
    nn_lsm  = 0             !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
                            !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
@@ -276 +273 @@
                            !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled
                            !  = 2  Redistribute a single flux over categories (coupled mode only)
+   nn_drag   = 0      !  formula to calculate momentum from the wind components 
+                           ! = 0 UKMO SHELF formulation 
+                           ! = 1 standard formulation with forced of coupled drag coefficient 
+                           ! = 2 standard formulation with constant drag coefficient 
+                           ! = 3 momentum calculated from core forcing fields
 /
 !-----------------------------------------------------------------------
@@ -298 +300 @@
    sn_emp      = 'emp'       ,        24         , 'emp'     , .false.      , .false., 'yearly'  , ''       , ''       , ''
 
-   cn_dir      = './'      !  root directory for the location of the flux files
+   cn_dir       = './'      !  root directory for the location of the flux files 
+   ln_shelf_flx = .false.   !  UKMO SHELF specific flux flag - read from file wind components instead of momentum  
+   ln_rel_wind  = .false.   !  UKMO SHELF - calculate momentum from wind speed relative to currents  
+   rn_wfac      = 1.0       !  UKMO SHELF - multiplicative factor for ocean/wind velocity
 /
 !-----------------------------------------------------------------------
@@ -383 +388 @@
    sn_rcv_mslp   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
    sn_rcv_phioc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
-   sn_rcv_sdrfx  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
-   sn_rcv_sdrfy  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_sdrft  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_wper   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_wfreq  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_wnum   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
-   sn_rcv_wstrf  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_tauoc  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''  
+   sn_rcv_tauw   =       'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
    sn_rcv_wdrag  =       'none'                 ,    'no'    ,     ''      ,         ''          ,   '' 
 !
@@ -936 +942 @@
    ln_zdfexp   = .false.   !  time-stepping: split-explicit (T) or implicit (F) time stepping
    nn_zdfexp   =    3            !  number of sub-timestep for ln_zdfexp=T
-   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010) (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)
 /
 !-----------------------------------------------------------------------
@@ -1000 +1005 @@
    rn_clim_galp  = 0.267   !  galperin limit
    ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case
-   rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux
+   rn_crban_default = 100.    !  Craig and Banner 1994 constant for wb tke flux
    rn_charn      = 70000.  !  Charnock constant for wb induced roughness length
    rn_hsro       =  0.02   !  Minimum surface roughness
@@ -1294 +1299 @@
 !              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
 !              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
+   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model 
+   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift 
+   ln_stcor    = .false.   !  Activate Stokes Coriolis term 
+   ln_tauoc    = .false.   !  Activate ocean stress modified by external wave induced stress 
+   ln_tauw     = .false.   !  Activate ocean stress components from wave model 
+   ln_phioc    = .false.   !  Activate wave to ocean energy 
+   ln_rough    = .false.   !  Wave roughness equals the significant wave height 
+   ln_zdfqiao  = .false.   !  Enhanced wave vertical mixing Qiao (2010) 
+   nn_sdrift   =  0   !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift 
+                           ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 
+                           ! = 1 Phillips:                      v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))] 
+   nn_wmix     =  0   !  type of wave breaking mixing 
+                           ! = 0 Craig and Banner formulation (original NEMO formulation)  
+                           ! = 1 Janssen formulation (no assumption of direct energy conversion)
    sn_cdg      =  'cdw_wave' ,        1          , 'drag_coeff' ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_usd      =  'sdw_wave' ,        1          , 'u_sd2d'     ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
@@ -1299 +1318 @@
    sn_hsw      =  'sdw_wave' ,        1          , 'hs'         ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
    sn_wmp      =  'sdw_wave' ,        1          , 'wmp'        ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_wfr      =  'sdw_wave' ,        1          , 'wave_freq'  ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''
    sn_wnum     =  'sdw_wave' ,        1          , 'wave_num'   ,     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
    sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_tauwx    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_tauwy    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   sn_phioc    =  'sdw_wave' ,        1          , 'wave_energy',     .true.   , .false. , 'daily'   ,  ''      , ''       , ''  
+   cn_dir      = './'  !  root directory for the location of wave forcing files
 !  
-   cn_dir  = './'  !  root directory for the location of drag coefficient files
 /
 !-----------------------------------------------------------------------

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

@@ -1094 +1094 @@
 
       ! Write all fields on T grid
+      IF( ln_wave .AND. ln_sdw ) THEN  
+         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity  
+         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity  
+         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity  
+      ENDIF
       CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
       CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
@@ -1111 +1116 @@
       END IF
        
-      IF( ln_wave .AND. ln_sdw ) THEN  
-         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity  
-         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity  
-         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity  
-      ENDIF  
-
       ! 3. Close the file
       ! -----------------

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -470 +470 @@
                 &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
       ENDIF
-      !  
-      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary  
-         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)  
-      ENDIF  
-      !  
 #if defined key_asminc
       !                                         ! Include the IAU weighted SSH increment
@@ -488 +483 @@
       !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations 
+      !  Phase 2 : Integration of the barotropic equations
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
-      !                                             ! ==================== !  
+      !                                             ! ==================== !
       ! Initialize barotropic variables:      
       IF( ll_init )THEN
@@ -502 +497 @@
          ub_e   (:,:) = 0._wp
          vb_e   (:,:) = 0._wp
+      ENDIF
+      ! 
+      ! Moved here to allow merging with other branch  
+      IF( ln_sdw ) THEN                         ! Stokes drift divergence added if necessary  
+         zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)  
       ENDIF
       !

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

@@ -182 +182 @@
          ELSE
                              CALL vor_een( kt, ntot, un, vn, ua, va )    ! total vorticity 
-            IF( ln_stcor )   CALL vor_ene( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
+            IF( ln_stcor )   CALL vor_een( kt, ncor, usd, vsd, ua, va )  ! add the Stokes-Coriolis trend
          ENDIF
          !

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

@@ -42 +42 @@
    LOGICAL , PUBLIC ::   ln_cpl         !: ocean-atmosphere coupled formulation
    LOGICAL , PUBLIC ::   ln_mixcpl      !: ocean-atmosphere forced-coupled mixed formulation
+   LOGICAL , PUBLIC ::   ln_wavcpl      !: ocean-wave forced-coupled mixed formulation 
+   LOGICAL , PUBLIC ::   ll_purecpl     !: ocean-atmosphere or ocean-wave pure coupled formulation
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
@@ -64 +66 @@
    LOGICAL , PUBLIC ::   ln_wave        !: true if some coupling with wave model
    LOGICAL , PUBLIC ::   ln_cdgw        !: true if neutral drag coefficient from wave model
-   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
+   LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d Stokes drift from wave model 
+   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
    LOGICAL , PUBLIC ::   ln_tauoc       !: true if normalized stress from wave is used  
-   LOGICAL , PUBLIC ::   ln_stcor       !: true if Stokes-Coriolis term is used
+   LOGICAL , PUBLIC ::   ln_tauw        !: true if ocean stress components from wave is used  
+   LOGICAL , PUBLIC ::   ln_phioc       !: true if wave energy to ocean is used  
+   LOGICAL , PUBLIC ::   ln_rough       !: true if wave roughness equals significant wave height 
+   LOGICAL , PUBLIC ::   ln_zdfqiao     !: Enhanced wave vertical mixing Qiao(2010) formulation flag 
+   INTEGER , PUBLIC ::   nn_drag        ! type of formula to calculate wind stress from wind components 
+   INTEGER , PUBLIC ::   nn_wmix        ! type of wave breaking mixing
    !
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
@@ -93 +101 @@
    INTEGER , PUBLIC, PARAMETER ::   jp_iam_sas  = 2      !: Multi executable configuration - SAS component
                                                          !  (internal OASIS coupling)
+   !!---------------------------------------------------------------------- 
+   !!           wind stress definition 
+   !!---------------------------------------------------------------------- 
+   INTEGER, PUBLIC, PARAMETER ::   jp_ukmo  = 0        ! UKMO SHELF formulation 
+   INTEGER, PUBLIC, PARAMETER ::   jp_std   = 1        ! standard formulation with forced or coupled drag coefficient  
+   INTEGER, PUBLIC, PARAMETER ::   jp_const = 2        ! standard formulation with constant drag coefficient  
+   INTEGER, PUBLIC, PARAMETER ::   jp_mcore = 3        ! momentum calculated from core forcing fields  
+
+   !!---------------------------------------------------------------------- 
+   !!           Wave mixing vertical parameterization 
+   !!---------------------------------------------------------------------- 
+   INTEGER, PUBLIC, PARAMETER ::   jp_craigbanner = 0  ! Craig and Banner formulation (original NEMO formulation - 
+                                                       ! direct conversion of mechanical to turbulent energy) 
+   INTEGER, PUBLIC, PARAMETER ::   jp_janssen     = 1  ! Janssen formulation - no assumption on direct energy conversion
    !!----------------------------------------------------------------------
    !!              Ocean Surface Boundary Condition fields

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90

@@ -26 +26 @@
    
    !                                !!* namsbc_apr namelist (Atmospheric PRessure) *
+   LOGICAL, PUBLIC ::   cpl_mslp = .FALSE. ! Is the pressure read from coupling?
    LOGICAL, PUBLIC ::   ln_apr_obc   !: inverse barometer added to OBC ssh data 
    LOGICAL, PUBLIC ::   ln_ref_apr   !: ref. pressure: global mean Patm (F) or a constant (F)
-   REAL(wp)        ::   rn_pref      !  reference atmospheric pressure   [N/m2]
+   REAL(wp),PUBLIC ::   rn_pref      !  reference atmospheric pressure   [N/m2]
 
    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   ssh_ib    ! Inverse barometer now    sea surface height   [m]
@@ -34 +35 @@
    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   apr       ! atmospheric pressure at kt                 [N/m2]
    
-   REAL(wp) ::   tarea                ! whole domain mean masked ocean surface
-   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
+   REAL(wp), PUBLIC ::   tarea                ! whole domain mean masked ocean surface
+   REAL(wp), PUBLIC ::   r1_grau              ! = 1.e0 / (grav * rau0)
    
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_apr   ! structure of input fields (file informations, fields read)
@@ -85 +86 @@
          IF(lwm) WRITE ( numond, namsbc_apr )
          !
-         ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst
-         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' )
-         !
-         CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' )
-                                ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   )
-         IF( sn_apr%ln_tint )   ALLOCATE( sf_apr(1)%fdta(jpi,jpj,1,2) )
+         IF( .NOT. cpl_mslp ) THEN 
+            ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst 
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' ) 
+            ! 
+            CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' ) 
+                                   ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   ) 
+            IF( sn_apr%ln_tint )   ALLOCATE( sf_apr(1)%fdta(jpi,jpj,1,2) ) 
+         ENDIF
                                 ALLOCATE( ssh_ib(jpi,jpj) , ssh_ibb(jpi,jpj) )
                                 ALLOCATE( apr (jpi,jpj) )
@@ -96 +99 @@
          IF(lwp) THEN                                 !* control print
             WRITE(numout,*)
-            WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric PRessure as extrenal forcing'
+            IF( cpl_mslp ) THEN 
+               WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric Pressure as extrenal coupling' 
+            ELSE 
+               WRITE(numout,*) '   Namelist namsbc_apr : Atmospheric Pressure as extrenal forcing' 
+            ENDIF
             WRITE(numout,*) '      ref. pressure: global mean Patm (T) or a constant (F)  ln_ref_apr = ', ln_ref_apr
          ENDIF
@@ -119 +126 @@
       ENDIF
 
-      !                                         ! ========================== !
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    At each sbc time-step   !
-         !                                      ! ===========+++============ !
-         !
-         IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields
-         !
-         CALL fld_read( kt, nn_fsbc, sf_apr )               !* input Patm provided at kt + nn_fsbc/2
-         !
-         !                                                  !* update the reference atmospheric pressure (if necessary)
-         IF( ln_ref_apr )   rn_pref = glob_sum( sf_apr(1)%fnow(:,:,1) * e1e2t(:,:) ) / tarea
-         !
-         !                                                  !* Patm related forcing at kt
-         ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)
-         apr   (:,:) =     sf_apr(1)%fnow(:,:,1)                        ! atmospheric pressure
-         !
-         CALL iom_put( "ssh_ib", ssh_ib )                   !* output the inverse barometer ssh
-      ENDIF
-
-      !                                         ! ---------------------------------------- !
-      IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    !
-         !                                      ! ---------------------------------------- !
-         !                                            !* Restart: read in restart file
-         IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN 
-            IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file'
-            CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh
+      IF( .NOT. cpl_mslp ) THEN                    ! ========================== ! 
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    At each sbc time-step   ! 
+            !                                      ! ===========+++============ !
             !
-         ELSE                                         !* no restart: set from nit000 values
-            IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb set to nit000 values'
-            ssh_ibb(:,:) = ssh_ib(:,:)
+            IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields 
+            ! 
+            CALL fld_read( kt, nn_fsbc, sf_apr )               !* input Patm provided at kt + nn_fsbc/2 
+            ! 
+            !                                                  !* update the reference atmospheric pressure (if necessary) 
+            IF( ln_ref_apr )   rn_pref = glob_sum( sf_apr(1)%fnow(:,:,1) * e1e2t(:,:) ) / tarea 
+            ! 
+            !                                                  !* Patm related forcing at kt 
+            ssh_ib(:,:) = - ( sf_apr(1)%fnow(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer) 
+            apr   (:,:) =     sf_apr(1)%fnow(:,:,1)                        ! atmospheric pressure 
+            ! 
+            CALL iom_put( "ssh_ib", ssh_ib )                   !* output the inverse barometer ssh
          ENDIF
-      ENDIF
-      !                                         ! ---------------------------------------- !
-      IF( lrst_oce ) THEN                       !      Write in the ocean restart file     !
-         !                                      ! ---------------------------------------- !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'sbc_apr : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp
-         IF(lwp) WRITE(numout,*) '~~~~'
-         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib )
+ 
+         !                                         ! ---------------------------------------- ! 
+         IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    ! 
+            !                                      ! ---------------------------------------- ! 
+            !                                            !* Restart: read in restart file 
+            IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN  
+               IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file' 
+               CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh 
+               ! 
+            ELSE                                         !* no restart: set from nit000 values 
+               IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb set to nit000 values' 
+               ssh_ibb(:,:) = ssh_ib(:,:) 
+            ENDIF 
+         ENDIF 
+         !                                         ! ---------------------------------------- ! 
+         IF( lrst_oce ) THEN                       !      Write in the ocean restart file     ! 
+            !                                      ! ---------------------------------------- ! 
+            IF(lwp) WRITE(numout,*) 
+            IF(lwp) WRITE(numout,*) 'sbc_apr : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp 
+            IF(lwp) WRITE(numout,*) '~~~~' 
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib ) 
+         ENDIF
       ENDIF
       !

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -319 +319 @@
          &               Cd, Ch, Ce, zt_zu, zq_zu )
     
+      IF ( ln_cdgw .AND. nn_drag==jp_std ) Cd(:,:) = cdn_wave(:,:)
       ! ... tau module, i and j component
       DO jj = 1, jpj
@@ -737 +738 @@
 
       !! Neutral coefficients at 10m:
-      IF( ln_wave .AND. ln_cdgw ) THEN      ! wave drag case
+      IF( ln_cdgw .AND. nn_drag==jp_mcore ) THEN      ! wave drag case
          cdn_wave(:,:) = cdn_wave(:,:) + rsmall * ( 1._wp - tmask(:,:,1) )
          ztmp0   (:,:) = cdn_wave(:,:)
@@ -783 +784 @@
          END IF
        
-         IF( ln_wave .AND. ln_cdgw ) THEN      ! surface wave case
+         IF( ln_cdgw .AND. nn_drag==jp_mcore ) THEN      ! surface wave case
             sqrt_Cd = vkarmn / ( vkarmn / sqrt_Cd_n10 - zpsi_m_u ) 
             Cd      = sqrt_Cd * sqrt_Cd

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -116 +116 @@
    INTEGER, PARAMETER ::   jpr_wper   = 48            ! Mean wave period  
    INTEGER, PARAMETER ::   jpr_wnum   = 49            ! Mean wavenumber  
-   INTEGER, PARAMETER ::   jpr_wstrf  = 50            ! Stress fraction adsorbed by waves  
+   INTEGER, PARAMETER ::   jpr_tauoc  = 50            ! Stress fraction adsorbed by waves
    INTEGER, PARAMETER ::   jpr_wdrag  = 51            ! Neutral surface drag coefficient  
-   INTEGER, PARAMETER ::   jprcv      = 51            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_wfreq  = 52            ! Wave peak frequency  
+   INTEGER, PARAMETER ::   jpr_tauwx  = 53            ! x component of the ocean stress from waves 
+   INTEGER, PARAMETER ::   jpr_tauwy  = 54            ! y component of the ocean stress from waves 
+   INTEGER, PARAMETER ::   jprcv      = 54            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -170 +173 @@
    TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev   
    ! Received from waves   
-   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrfx,sn_rcv_sdrfy,sn_rcv_wper,sn_rcv_wnum,sn_rcv_wstrf,sn_rcv_wdrag
+   TYPE(FLD_C) ::   sn_rcv_hsig,sn_rcv_phioc,sn_rcv_sdrft,sn_rcv_wper, & 
+                    sn_rcv_wfreq,sn_rcv_wnum,sn_rcv_tauoc,sn_rcv_tauw, & 
+                    sn_rcv_wdrag
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -183 +188 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
     
-   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]   
-   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0
-
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
 
@@ -247 +249 @@
          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,      &   
          &                  sn_snd_ifrac, sn_snd_crtw , sn_snd_wlev  , sn_rcv_hsig  , sn_rcv_phioc ,   &   
-         &                  sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wstrf ,   &  
-         &                  sn_rcv_wdrag, sn_rcv_qns  , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   ,   &  
-         &                  sn_rcv_iceflx,sn_rcv_co2  , nn_cplmodel  , ln_usecplmask, sn_rcv_mslp
+         &                  sn_rcv_sdrft, sn_rcv_wper  , sn_rcv_wnum  , sn_rcv_wfreq, sn_rcv_tauoc,    & 
+         &                  sn_rcv_wdrag, sn_rcv_qns   , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal ,     & 
+         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp , sn_rcv_tauw ,                  & 
+         &                  nn_cplmodel, ln_usecplmask
       !!---------------------------------------------------------------------
       !
@@ -290 +293 @@
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
+         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')'
          WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')'   
          WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')'   
-         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')'   
-         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')'   
+         WRITE(numout,*)'      Surface Stokes drift u,v        = ', TRIM(sn_rcv_sdrft%cldes ), ' (', TRIM(sn_rcv_sdrft%clcat ), ')'
          WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')'   
          WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')'   
-         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_wstrf%cldes ), ' (', TRIM(sn_rcv_wstrf%clcat ), ')'   
+         WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'   
+         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauoc%cldes ), ' (', TRIM(sn_rcv_tauoc%clcat ), ')'   
+         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
          WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
@@ -309 +314 @@
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
          WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')'   
-         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')'   
          WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')'   
          WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref   
@@ -525 +529 @@
       !                                                      ! Mean Sea Level Pressure   !   
       !                                                      ! ------------------------- !   
-      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE.   
+      srcv(jpr_mslp)%clname = 'O_MSLP' 
+      IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' ) THEN 
+         srcv(jpr_mslp)%laction = .TRUE. 
+         cpl_mslp = .TRUE. 
+      ENDIF
 
       !                                                      ! ------------------------- !
@@ -554 +562 @@
       ENDIF  
       srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction  
-      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN  
+      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction  
+      IF( TRIM(sn_rcv_sdrft%cldes ) == 'coupled' )  THEN
          srcv(jpr_sdrftx)%laction = .TRUE.  
-         cpl_sdrftx = .TRUE.  
-      ENDIF  
-      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction  
-      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN  
-         srcv(jpr_sdrfty)%laction = .TRUE.  
-         cpl_sdrfty = .TRUE.  
+         cpl_sdrft = .TRUE.
       ENDIF  
       srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period  
@@ -568 +572 @@
          cpl_wper = .TRUE.  
       ENDIF  
+      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency  
+      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN  
+         srcv(jpr_wfreq)%laction = .TRUE.  
+         cpl_wfreq = .TRUE.  
+      ENDIF
       srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number  
       IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN  
@@ -573 +582 @@
          cpl_wnum = .TRUE.  
       ENDIF  
-      srcv(jpr_wstrf)%clname = 'O_WStrf'     ! stress fraction adsorbed by the wave  
-      IF( TRIM(sn_rcv_wstrf%cldes ) == 'coupled' )  THEN  
-         srcv(jpr_wstrf)%laction = .TRUE.  
-         cpl_wstrf = .TRUE.  
+      srcv(jpr_tauoc)%clname = 'O_TauOce'     ! stress fraction adsorbed by the wave  
+      IF( TRIM(sn_rcv_tauoc%cldes ) == 'coupled' )  THEN  
+         srcv(jpr_tauoc)%laction = .TRUE.  
+         cpl_tauoc = .TRUE.  
+      ENDIF  
+      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction 
+      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction 
+      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN  
+         srcv(jpr_tauwx)%laction = .TRUE.  
+         srcv(jpr_tauwy)%laction = .TRUE.  
+         cpl_tauw = .TRUE.
       ENDIF  
       srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient  
@@ -583 +599 @@
          cpl_wdrag = .TRUE.  
       ENDIF  
+      ! 
+      IF( srcv(jpr_tauoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) & 
+            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', & 
+                                     '(sn_rcv_tauoc=coupled and sn_rcv_tauw=coupled)' ) 
+      !
       !   
       !                                                      ! ------------------------------- !
@@ -964 +985 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      USE zdf_oce,  ONLY : ln_zdfqiao
+      USE sbcflx ,  ONLY : ln_shelf_flx 
+      USE sbcssm ,  ONLY : sbc_ssm_cpl 
+      USE lib_fortran     ! distributed memory computing library
 
       INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
@@ -1071 +1094 @@
       ELSE                                                   !   No dynamical coupling   !
          !                                                   ! ========================= !
+         ! it is possible that the momentum is calculated from the winds (ln_shelf_flx) and a coupled drag coefficient 
+         IF( srcv(jpr_wdrag)%laction .AND. ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) THEN 
+            DO jj = 1, jpjm1 
+               DO ji = 1, jpim1 
+                  ! here utau and vtau should contain the wind components as read from the forcing files, 
+                  ! and the wind module should already be properly calculated 
+                  frcv(jpr_otx1)%z3(ji,jj,1) = zrhoa * 0.5 * ( frcv(jpr_wdrag)%z3(ji,jj,1) + frcv(jpr_wdrag)%z3(ji+1,jj,1) ) * & 
+                                                                         utau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
+                  frcv(jpr_oty1)%z3(ji,jj,1) = zrhoa * 0.5 * ( frcv(jpr_wdrag)%z3(ji,jj,1) + frcv(jpr_wdrag)%z3(ji,jj+1,1) ) * & 
+                                                                         vtau(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
+                  utau(ji,jj) = frcv(jpr_otx1)%z3(ji,jj,1) 
+                  vtau(ji,jj) = frcv(jpr_oty1)%z3(ji,jj,1) 
+               END DO 
+            END DO 
+            CALL lbc_lnk_multi( frcv(jpr_otx1)%z3(:,:,1), 'U', -1. , frcv(jpr_oty1)%z3(:,:,1), 'V', -1. , & 
+                                                             utau(:,:), 'U', -1. , vtau(:,:), 'V',  -1. ) 
+            llnewtx = .TRUE. 
+         ELSE
          frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero 
          frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
          llnewtx = .TRUE.
+         ENDIF
          !
       ENDIF
@@ -1093 +1135 @@
             END DO
             CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
+            IF( .NOT. srcv(jpr_otx1)%laction .AND. srcv(jpr_wdrag)%laction .AND. & 
+                                ln_shelf_flx .AND. ln_cdgw .AND. nn_drag == jp_std ) & 
+               taum(:,:) = frcv(jpr_taum)%z3(:,:,1
             llnewtau = .TRUE.
          ELSE
@@ -1107 +1152 @@
       !                                                      ! ========================= !
       !                                                      !      10 m wind speed      !   (wndm)
+      !                                                      !   include wave drag coef  !   (wndm)
       !                                                      ! ========================= !
       !
@@ -1117 +1163 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi 
+                  IF( ln_shelf_flx ) THEN   ! the 10 wind module is properly calculated before if ln_shelf_flx 
+                     frcv(jpr_w10m)%z3(ji,jj,1) = wndm(ji,jj) 
+                  ELSE
                   frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
+                  ENDIF
                END DO
             END DO
@@ -1127 +1177 @@
       IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
          !
+         ! if ln_wavcpl, the fields already contain the right information from forcing even if not ln_mixcpl
          IF( ln_mixcpl ) THEN
-            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
-            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
-            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
-            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
-         ELSE
+            IF( srcv(jpr_otx1)%laction ) THEN 
+               utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:) 
+               vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:) 
+            ENDIF 
+            IF( srcv(jpr_taum)%laction )   & 
+               taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:) 
+            IF( srcv(jpr_w10m)%laction )   & 
+               wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN
             utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
             vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
@@ -1156 +1211 @@
           IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields   
        
-          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization   
-          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)   
+          !                                                  !* update the reference atmospheric pressure (if necessary) 
+          IF( ln_ref_apr )  rn_pref = glob_sum( frcv(jpr_mslp)%z3(:,:,1) * e1e2t(:,:) ) / tarea 
+        
+          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rn_pref ) * r1_grau    ! equivalent ssh (inverse barometer)   
           apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1) !atmospheric pressure   
+          ! 
+          CALL iom_put( "ssh_ib", ssh_ib )                                    !* output the inverse barometer ssh
        
-          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)   
+          !                                         ! ---------------------------------------- ! 
+          IF( kt == nit000 ) THEN                   !   set the forcing field at nit000 - 1    ! 
+             !                                      ! ---------------------------------------- ! 
+             !* Restart: read in restart file 
+             IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN 
+                IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb read in the restart file' 
+                CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb )   ! before inv. barometer ssh 
+             ELSE                                         !* no restart: set from nit000 values 
+                IF(lwp) WRITE(numout,*) 'sbc_cpl:   ssh_ibb set to nit000 values' 
+                ssh_ibb(:,:) = ssh_ib(:,:) 
+             ENDIF 
+          ENDIF 
+          !                                         ! ---------------------------------------- ! 
+          IF( lrst_oce ) THEN                       !      Write in the ocean restart file     ! 
+             !                                      ! ---------------------------------------- ! 
+             IF(lwp) WRITE(numout,*) 
+             IF(lwp) WRITE(numout,*) 'sbc_cpl : ssh_ib written in ocean restart file at it= ', kt,' date= ', ndastp 
+             IF(lwp) WRITE(numout,*) '~~~~' 
+             CALL iom_rstput( kt, nitrst, numrow, 'ssh_ibb' , ssh_ib ) 
+          ENDIF 
+        
+          ! Update mean ssh 
+          CALL sbc_ssm_cpl( kt )
       END IF   
       !  
       IF( ln_sdw ) THEN  ! Stokes Drift correction activated  
       !                                                      ! ========================= !   
-      !                                                      !       Stokes drift u      !  
+      !                                                      !       Stokes drift u,v    !  
       !                                                      ! ========================= !   
-         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)  
-      !  
-      !                                                      ! ========================= !   
-      !                                                      !       Stokes drift v      !  
-      !                                                      ! ========================= !   
-         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)  
+      IF( srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction ) THEN 
+                                     ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)  
+                                     vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)  
+      ENDIF
       !  
       !                                                      ! ========================= !   
@@ -1185 +1264 @@
       !  
       !                                                      ! ========================= !   
+      !                                                      !    Wave peak frequency    !  
+      !                                                      ! ========================= !   
+         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)  
+      ! 
+      !                                                      ! ========================= !
       !                                                      !    Vertical mixing Qiao   !  
       !                                                      ! ========================= !   
@@ -1190 +1274 @@
        
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode  
-         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &  
-                                                                    .OR. srcv(jpr_hsig)%laction ) &  
+         IF( (srcv(jpr_sdrftx)%laction .AND. srcv(jpr_sdrfty)%laction) .OR. srcv(jpr_wper)%laction &  
+                                        .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) &
             CALL sbc_stokes()  
       ENDIF  
@@ -1197 +1281 @@
       !                                                      ! Stress adsorbed by waves  !  
       !                                                      ! ========================= !   
-      IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1)  
+      IF( srcv(jpr_tauoc)%laction .AND. ln_tauoc ) THEN 
+         tauoc_wave(:,:) = frcv(jpr_tauoc)%z3(:,:,1) 
+         ! cap the value of tauoc 
+         WHERE(tauoc_wave <   0.0 ) tauoc_wave = 1.0 
+         WHERE(tauoc_wave > 100.0 ) tauoc_wave = 1.0 
+      ENDIF 
+      !                                                      ! ========================= !   
+      !                                                      ! Stress component by waves !  
+      !                                                      ! ========================= !   
+      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN 
+         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1) 
+         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1) 
+         ! cap the value of tauoc 
+         WHERE(tauw_x < -100.0 ) tauw_x = 0.0 
+         WHERE(tauw_x >  100.0 ) tauw_x = 0.0 
+         WHERE(tauw_y < -100.0 ) tauw_y = 0.0 
+         WHERE(tauw_y >  100.0 ) tauw_y = 0.0 
+      ENDIF
        
       !                                                      ! ========================= !   
-      !                                                      ! Wave drag coefficient   !  
+      !                                                      !   Wave to ocean energy    !
       !                                                      ! ========================= !   
-      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
+      IF( srcv(jpr_phioc)%laction .AND. ln_phioc ) THEN 
+         rn_crban(:,:) = 29.0 * frcv(jpr_phioc)%z3(:,:,1) 
+         WHERE( rn_crban <    0.0 ) rn_crban = 0.0 
+         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0 
+      ENDIF 
 
       !  Fields received by SAS when OASIS coupling
@@ -1274 +1379 @@
                CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
             END SELECT
-         ELSE
+         ELSE IF( ll_purecpl ) THEN
             zemp(:,:) = 0._wp
          ENDIF
@@ -1282 +1387 @@
          IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
          
-         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
-         ELSE                   ;   emp(:,:) =                              zemp(:,:)
+         IF( ln_mixcpl .AND. ( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction )) THEN 
+                                         emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN  ;   emp(:,:) = zemp(:,:)
          ENDIF
          !
@@ -1298 +1404 @@
             ENDIF
          ENDIF
-         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
-         ELSE                   ;   qns(:,:) =                              zqns(:,:)
+         IF( ln_mixcpl .AND. ( srcv(jpr_qnsoce)%laction .OR. srcv(jpr_qnsmix)%laction )) THEN 
+                                          qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN   ;   qns(:,:) = zqns(:,:)
          ENDIF
 
@@ -1308 +1415 @@
          ENDIF
          IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
-         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
-         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
+         IF( ln_mixcpl .AND. ( srcv(jpr_qsroce)%laction .OR. srcv(jpr_qsrmix)%laction )) THEN 
+                                          qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:) 
+         ELSE IF( ll_purecpl ) THEN   ;   qsr(:,:) = zqsr(:,:)
          ENDIF
          !

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

@@ -20 +20 @@
    USE iom             ! IOM library
    USE in_out_manager  ! I/O manager
+   USE sbcwave         ! wave physics
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -87 +88 @@
       REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
+      REAL(wp) ::   totwind               ! UKMO SHELF: Module of wind speed
       REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
       REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
@@ -115 +117 @@
          IF(lwm) WRITE ( numond, namsbc_flx ) 
          !
+         IF(lwp) THEN                        ! Namelist print 
+            WRITE(numout,*) 
+            WRITE(numout,*) 'sbc_flx : Flux forcing' 
+            WRITE(numout,*) '~~~~~~~~~~~' 
+            WRITE(numout,*) '       Namelist namsbc_flx : shelf seas configuration (force with winds instead of momentum)' 
+            WRITE(numout,*) '          shelf seas configuration    ln_shelf_flx    = ', ln_shelf_flx 
+            WRITE(numout,*) '          relative wind speed         ln_rel_wind     = ', ln_rel_wind 
+            WRITE(numout,*) '          wind multiplication factor  rn_wfac         = ', rn_wfac 
+         ENDIF
          !                                         ! check: do we plan to use ln_dm2dc with non-daily forcing?
          IF( ln_dm2dc .AND. sn_qsr%nfreqh /= 24 )   &
@@ -210 +221 @@
             END DO
          END DO
+         !                                                        ! add modification due to drag coefficient read from wave forcing 
+         !                                                        ! this code is inefficient but put here to allow merging with another UKMO branch 
+         IF( ln_shelf_flx ) THEN 
+            ! calculate first the wind module, as it will be used later 
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1    ! vect. opt. 
+                  ztx = zwnd_i(ji-1,jj  ) + zwnd_i(ji,jj) 
+                  zty = zwnd_j(ji  ,jj-1) + zwnd_j(ji,jj) 
+                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty ) 
+               END DO 
+            END DO 
+            CALL lbc_lnk( wndm(:,:), 'T', 1. ) 
+        
+            IF( ln_cdgw .AND. nn_drag == jp_std ) THEN 
+               IF( cpl_wdrag ) THEN  
+                  ! reset utau and vtau to the wind components: the momentum will 
+                  ! be calculated from the coupled value of the drag coefficient 
+                  DO jj = 1, jpj 
+                     DO ji = 1, jpi 
+                        utau(ji,jj) = zwnd_i(ji,jj) 
+                        vtau(ji,jj) = zwnd_j(ji,jj) 
+                     END DO 
+                  END DO 
+               ELSE 
+                  DO jj = 1, jpjm1 
+                     DO ji = 1, jpim1 
+                        utau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * & 
+                                                                        0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
+                        vtau(ji,jj) = zrhoa * 0.5 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * & 
+                                                                        0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
+                     END DO 
+                  END DO 
+                  CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
+               ENDIF 
+            ELSE IF( nn_drag == jp_const ) THEN 
+               DO jj = 1, jpjm1 
+                  DO ji = 1, jpim1 
+                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji+1,jj) ) 
+                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * 0.5 * ( wndm(ji,jj) + wndm(ji,jj+1) ) 
+                  END DO 
+               END DO 
+               CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
+            ENDIF 
+         ENDIF
          !                                                        ! add to qns the heat due to e-p
          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
@@ -230 +285 @@
                zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
                taum(ji,jj) = zmod
+               IF( .NOT. ln_shelf_flx ) THEN
                wndm(ji,jj) = SQRT( zmod * zcoef )
+               ENDIF
             END DO
          END DO

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

@@ -88 +88 @@
       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
-         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
-         &             ln_tauoc  , ln_stcor  , nn_lsm, nn_limflx , nn_components, ln_cpl
+         &             ln_ssr    , nn_isf    , nn_fwb, ln_wave, nn_lsm     , nn_limflx,   & 
+                       nn_components, ln_cpl , ln_wavcpl, nn_drag
       INTEGER  ::   ios
       INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
-      LOGICAL  ::   ll_purecpl
       !!----------------------------------------------------------------------
 
@@ -131 +130 @@
          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
          WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
-         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
+         WRITE(numout,*) '              forced-coupled atm mixed formulation       ln_mixcpl   = ', ln_mixcpl 
+         WRITE(numout,*) '              forced-coupled wav mixed formulation       ln_wavcpl   = ', ln_wavcpl
          WRITE(numout,*) '              wave physics                               ln_wave     = ', ln_wave  
-         WRITE(numout,*) '                 Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw  
-         WRITE(numout,*) '                 wave modified ocean stress              ln_tauoc    = ', ln_tauoc  
-         WRITE(numout,*) '                 Stokes coriolis term                    ln_stcor    = ', ln_stcor  
-         WRITE(numout,*) '                 neutral drag coefficient (CORE, MFS)    ln_cdgw     = ', ln_cdgw
          WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
          WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
          WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
+         WRITE(numout,*) '              momentum formulation                       nn_drag     = ', nn_drag
          WRITE(numout,*) '           Misc. options of sbc : '
          WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
@@ -169 +166 @@
       IF ( nn_components == jp_iam_opa .AND. ln_cpl )   &
          &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
-      IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
-         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      IF ( nn_components == jp_iam_opa .AND. ( ln_mixcpl .OR. ln_wavcpl) )   & 
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl or ln_wavcpl = T in OPA' )
       IF ( ln_cpl .AND. .NOT. lk_oasis )    &
          &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
-      IF( ln_mixcpl .AND. .NOT. lk_oasis )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' )
-      IF( ln_mixcpl .AND. .NOT. ln_cpl )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' )
-      IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    &
-         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' )
+      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. lk_oasis )    & 
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires the cpp key key_oasis3' ) 
+      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. .NOT. ln_cpl )    & 
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) requires ln_cpl = T' ) 
+      IF( ( ln_mixcpl .OR. ln_wavcpl ) .AND. nn_components /= jp_iam_nemo )    & 
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl or ln_wavcpl) is not yet working with sas-opa coupling via oasis' )
 
       !                              ! allocate sbc arrays
@@ -219 +216 @@
       IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   &
          &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
-      
-      IF ( ln_wave ) THEN
-         !Activated wave module but neither drag nor stokes drift activated  
-         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) )   THEN   
-             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_stcor=F')  
-         !drag coefficient read from wave model definable only with mfs bulk formulae and core
-         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
-             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
-         ELSEIF (ln_stcor .AND. .NOT. ln_sdw) THEN   
-             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
-         ENDIF
-      ELSE
-         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_stcor ) &   
-            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &  
-            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &  
-            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &  
-            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      &  
-            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
-      ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
-      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
+      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl .AND. .NOT. ln_wavcpl
       !
       icpt = 0
@@ -271 +249 @@
          IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
          IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
-         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation'
+         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed atm formulation' 
+         IF( ln_wavcpl          )   WRITE(numout,*) '              + forced-coupled mixed wav formulation'
          IF( nn_components/= jp_iam_nemo )  &
             &                       WRITE(numout,*) '              + OASIS coupled SAS'
@@ -340 +319 @@
       !!              - updte the ice fraction : fr_i
       !!----------------------------------------------------------------------
+      USE bdydta, ONLY: bdy_dta  
+      !
       INTEGER, INTENT(in) ::   kt       ! ocean time step
       !!---------------------------------------------------------------------
@@ -360 +341 @@
       !                                            ! ---------------------------------------- !
       !
-      IF ( .NOT. lk_bdy ) then
-         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
-      ENDIF
+
+      IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
                                                          ! (caution called before sbc_ssm)
       !
@@ -397 +377 @@
       END SELECT
 
-      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
-
-      IF ( ln_wave .AND. ln_tauoc) THEN                 ! Wave stress subctracted  
-            utau(:,:) = utau(:,:)*tauoc_wave(:,:)  
-            vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)  
-            taum(:,:) = taum(:,:)*tauoc_wave(:,:)  
-      !  
-            SELECT CASE( nsbc )  
-            CASE(  0,1,2,3,5,-1 )  ;  
-                IF(lwp .AND. kt == nit000 ) WRITE(numout,*) 'WARNING: You are subtracting the wave stress to the ocean. &  
-                        & If not requested select ln_tauoc=.false'  
-            END SELECT  
-      !  
-      END IF
+      IF( ln_mixcpl .OR. ln_wavcpl )  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
+      
+      IF( ln_wave .AND. (ln_tauoc .OR. ln_tauw) ) CALL sbc_stress( )    ! Wave stress update  
+      IF( lk_bdy )           CALL bdy_dta ( kt, time_offset=+1 )        ! update dynamic & tracer data at open boundaries 
+                                                                        ! (caution called after sbc_ssm[_cpl] and before ice)
 
       !                                            !==  Misc. Options  ==!

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

@@ -26 +26 @@
 
    PUBLIC   sbc_ssm         ! routine called by step.F90
+   PUBLIC   sbc_ssm_cpl     ! routine called by sbccpl.F90
    PUBLIC   sbc_ssm_init    ! routine called by sbcmod.F90
 
@@ -77 +78 @@
          sss_m(:,:) = zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = sshn(:,:)
+         IF( .NOT. cpl_mslp ) THEN 
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+            ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
+            ENDIF
          ENDIF
          !
@@ -99 +102 @@
             sss_m(:,:) = zcoef * zts(:,:,jp_sal)
             !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
-            ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:)
+            IF( .NOT. cpl_mslp ) THEN 
+               IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ) 
+               ELSE                    ;   ssh_m(:,:) = zcoef * sshn(:,:) 
+               ENDIF
             ENDIF
             !
@@ -113 +118 @@
             sst_m(:,:) = 0.e0
             sss_m(:,:) = 0.e0
-            ssh_m(:,:) = 0.e0
+            IF( .NOT. cpl_mslp ) ssh_m(:,:) = 0.e0
             IF( lk_vvl )   e3t_m(:,:) = 0.e0
             frq_m(:,:) = 0.e0
@@ -127 +132 @@
          sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
-         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
-         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:)
+         IF( .NOT. cpl_mslp ) THEN 
+            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+            ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
+            ENDIF
          ENDIF
          !
@@ -143 +150 @@
             ssu_m(:,:) = ssu_m(:,:) * zcoef           ! mean suface current  [m/s]
             ssv_m(:,:) = ssv_m(:,:) * zcoef           !
-            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
+            IF( .NOT. cpl_mslp ) ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
             IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m]
             frq_m(:,:) = frq_m(:,:) * zcoef   ! mean fraction of solar net radiation absorbed in the 1st T level [-]
@@ -161 +168 @@
             CALL iom_rstput( kt, nitrst, numrow, 'sst_m'  , sst_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  )
-            CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
+            IF( .NOT. cpl_mslp ) CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
             IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , e3t_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'frq_m'  , frq_m  )
@@ -174 +181 @@
          CALL iom_put( 'sst_m', sst_m )
          CALL iom_put( 'sss_m', sss_m )
-         CALL iom_put( 'ssh_m', ssh_m )
+         IF( .NOT. cpl_mslp ) CALL iom_put( 'ssh_m', ssh_m )
          IF( lk_vvl )   CALL iom_put( 'e3t_m', e3t_m )
          CALL iom_put( 'frq_m', frq_m )
@@ -180 +187 @@
       !
    END SUBROUTINE sbc_ssm
+
+   SUBROUTINE sbc_ssm_cpl( kt ) 
+      !!--------------------------------------------------------------------- 
+      !!                   ***  ROUTINE sbc_ssm_cpl  *** 
+      !!                      
+      !! ** Purpose :   provide ocean surface variable to sea-surface boundary 
+      !!                condition computation when pressure is read from coupling 
+      !!                 
+      !! ** Method  :   The inverse barometer ssh (i.e. ssh associated with Patm) 
+      !!                is added to ssh_m when ln_apr_dyn = T. Required for sea-ice dynamics. 
+      !!--------------------------------------------------------------------- 
+      INTEGER, INTENT(in) ::   kt   ! ocean time step 
+      ! 
+      REAL(wp) ::   zcoef       ! local scalar 
+      !!--------------------------------------------------------------------- 
+      ! 
+      IF( nn_fsbc == 1 ) THEN                             !      Instantaneous surface fields        ! 
+         !                                                ! ---------------------------------------- ! 
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+         ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
+         ENDIF 
+      ELSE 
+         !                                                ! ----------------------------------------------- ! 
+         IF( kt == nit000 .AND. .NOT. l_ssm_mean ) THEN   !   Initialisation: 1st time-step, no input means ! 
+            !                                             ! ----------------------------------------------- ! 
+            IF(lwp) WRITE(numout,*) 
+            IF(lwp) WRITE(numout,*) '~~~~~~~   mean ssh field initialised to instantaneous values' 
+            zcoef = REAL( nn_fsbc - 1, wp ) 
+            zcoef = REAL( nn_fsbc - 1, wp ) 
+            IF( ln_apr_dyn ) THEN    ;  ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ) 
+            ELSE                     ;  ssh_m(:,:) = zcoef * sshn(:,:) 
+            ENDIF 
+            !                                             ! ---------------------------------------- ! 
+         ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   ! 
+            !                                             ! ---------------------------------------- ! 
+            ssh_m(:,:) = 0.e0 
+         ENDIF 
+   
+         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
+         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
+         ENDIF 
+         !                                                ! ---------------------------------------- ! 
+         IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   ! 
+            !                                             ! ---------------------------------------- ! 
+            zcoef = 1. / REAL( nn_fsbc, wp ) 
+            ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH [m] 
+         ENDIF 
+         !                                                ! ---------------------------------------- ! 
+         IF( lrst_oce ) THEN                              !      Write in the ocean restart file     ! 
+            !                                             ! ---------------------------------------- ! 
+            IF(lwp) WRITE(numout,*) 
+            IF(lwp) WRITE(numout,*) 'sbc_ssm_cpl : ssh mean field written in ocean restart file ',   & 
+               &                    'at it= ', kt,' date= ', ndastp 
+            IF(lwp) WRITE(numout,*) '~~~~~~~' 
+            CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  ) 
+         ENDIF 
+      ENDIF 
+      ! 
+      IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   ! 
+         CALL iom_put( 'ssh_m', ssh_m ) 
+      ENDIF 
+      ! 
+   END SUBROUTINE sbc_ssm_cpl
 
    SUBROUTINE sbc_ssm_init

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90

@@ -18 +18 @@
    USE oce            ! ocean variables
    USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE zdf_oce,  ONLY : ln_zdfqiao
    USE bdy_oce        ! open boundary condition variables
    USE domvvl         ! domain: variable volume layers
@@ -33 +32 @@
 
    PUBLIC   sbc_stokes      ! routine called in sbccpl 
+   PUBLIC   sbc_stress      ! routine called in sbcmod
    PUBLIC   sbc_wave        ! routine called in sbcmod 
    PUBLIC   sbc_wave_init   ! routine called in sbcmod
@@ -39 +39 @@
    LOGICAL, PUBLIC ::   cpl_hsig   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_phioc  = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrftx = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_sdrfty = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_sdrft  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wper   = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_wfreq  = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wnum   = .FALSE.
-   LOGICAL, PUBLIC ::   cpl_wstrf  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauoc  = .FALSE.
+   LOGICAL, PUBLIC ::   cpl_tauw   = .FALSE.
    LOGICAL, PUBLIC ::   cpl_wdrag  = .FALSE.
+
+   INTEGER ::   nn_sdrift      ! type of parameterization to calculate vertical Stokes drift
+   INTEGER, PARAMETER ::   jp_breivik  = 0     ! Breivik 2015: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
+   INTEGER, PARAMETER ::   jp_phillips = 1     ! Phillips:     v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]
+   INTEGER, PARAMETER ::   jp_peakph   = 2     ! Phillips using the peak wave number read from wave model instead of the inverse depth scale
 
    INTEGER ::   jpfld    ! number of files to read for stokes drift
@@ -51 +57 @@
    INTEGER ::   jp_hsw   ! index of significant wave hight      (m)      at T-point
    INTEGER ::   jp_wmp   ! index of mean wave period            (s)      at T-point
+   INTEGER ::   jp_wfr   ! index of wave peak frequency         (s^-1)   at T-point
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_cd    ! structure of input fields (file informations, fields read) Drag Coefficient
@@ -56 +63 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_wn    ! structure of input fields (file informations, fields read) wave number for Qiao
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauoc ! structure of input fields (file informations, fields read) normalized wave stress into the ocean
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_tauw ! structure of input fields (file informations, fields read) ocean stress components from wave model
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::  sf_phioc ! structure of input fields (file informations, fields read) wave to ocean energy 
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   cdn_wave            !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   hsw, wmp, wnum      !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   wfreq               !: 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   rn_crban            !: Craig and Banner constant for surface breaking waves mixing
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauoc_wave          !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_x              !:
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tauw_y              !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d               !:
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd              !: barotropic stokes drift divergence
@@ -64 +77 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp.
 
-   !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -98 +110 @@
       CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd ) 
       !
-      zfac = 2.0_wp * rpi / 16.0_wp
-      DO jj = 1, jpj               ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) )
-         DO ji = 1, jpi
-            ! Stokes drift velocity estimated from Hs and Tmean
-            ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
-            ! Stokes surface speed
-            tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
-            ! Wavenumber scale
-            zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
-         END DO
-      END DO
-      DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
-         DO ji = 1, jpim1
-            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
-            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
-            !
-            zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
-            zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
-         END DO
-      END DO
+      ! select parameterization for the calculation of vertical Stokes drift
+      ! exp. wave number at t-point
+      IF( nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips ) THEN   ! (Eq. (19) in Breivick et al. (2014) )
+         zfac = 2.0_wp * rpi / 16.0_wp
+         DO jj = 1, jpj               
+            DO ji = 1, jpi
+               ! Stokes drift velocity estimated from Hs and Tmean
+               ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp )
+               ! Stokes surface speed
+               tsd2d(ji,jj) = SQRT( ut0sd(ji,jj)*ut0sd(ji,jj) + vt0sd(ji,jj)*vt0sd(ji,jj))
+               ! Wavenumber scale
+               zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp )
+            END DO
+         END DO
+         DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
+            DO ji = 1, jpim1
+               zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) )
+               zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) )
+               !
+               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+            END DO
+         END DO
+      ELSE IF( nn_sdrift==jp_peakph ) THEN    ! peak wave number calculated from the peak frequency received by the wave model
+         DO jj = 1, jpjm1              
+            DO ji = 1, jpim1
+               zk_u(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji+1,jj)*wfreq(ji+1,jj) ) / grav
+               zk_v(ji,jj) = 0.5_wp * ( wfreq(ji,jj)*wfreq(ji,jj) + wfreq(ji,jj+1)*wfreq(ji,jj+1) ) / grav
+               !
+               zu0_sd(ji,jj) = 0.5_wp * ( ut0sd(ji,jj) + ut0sd(ji+1,jj) )
+               zv0_sd(ji,jj) = 0.5_wp * ( vt0sd(ji,jj) + vt0sd(ji,jj+1) )
+            END DO
+         END DO
+      ENDIF
       !
       !                       !==  horizontal Stokes Drift 3D velocity  ==!
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
-               zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
-               !                          
-               zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
-               zkh_v = zk_v(ji,jj) * zdep_v
-               !                                ! Depth attenuation
-               zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
-               zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
-               !
-               usd(ji,jj,jk) = zda_u * zk_u(ji,jj) * umask(ji,jj,jk)
-               vsd(ji,jj,jk) = zda_v * zk_v(ji,jj) * vmask(ji,jj,jk)
+      IF( nn_sdrift==jp_breivik ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
+                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
+                  !                          
+                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
+                  zkh_v = zk_v(ji,jj) * zdep_v
+                  !                                ! Depth attenuation
+                  zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
+                  zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
+                  !
+                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+               END DO
             END DO
          END DO
-      END DO
+      ELSE IF( nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
+                  zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
+                  !                          
+                  zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
+                  zkh_v = zk_v(ji,jj) * zdep_v
+                  !                                ! Depth attenuation
+                  zda_u = EXP( -2.0_wp*zkh_u ) - SQRT(2.0_wp*rpi*zkh_u) * ERFC(SQRT(2.0_wp*zkh_u))
+                  zda_v = EXP( -2.0_wp*zkh_v ) - SQRT(2.0_wp*rpi*zkh_v) * ERFC(SQRT(2.0_wp*zkh_v))
+                  !
+                  usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
+                  vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
       CALL lbc_lnk( usd(:,:,:), 'U', vsd(:,:,:), 'V', -1. )
       !
@@ -190 +237 @@
 
 
+   SUBROUTINE sbc_stress( )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_stress  ***
+      !!
+      !! ** Purpose :   Updates the ocean momentum modified by waves
+      !!
+      !! ** Method  : - Calculate u,v components of stress depending on stress
+      !!                model 
+      !!              - Calculate the stress module
+      !!              - The wind module is not modified by waves 
+      !! ** action  
+      !!---------------------------------------------------------------------
+      INTEGER  ::   jj, ji   ! dummy loop argument
+      !
+      IF( ln_tauoc ) THEN
+         utau(:,:) = utau(:,:)*tauoc_wave(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wave(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wave(:,:)
+      ENDIF
+      !
+      IF( ln_tauw ) THEN
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! Stress components at u- & v-points
+               utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) )
+               vtau(ji,jj) = 0.5_wp * ( tauw_y(ji,jj) + tauw_y(ji,jj+1) )
+               !
+               ! Stress module at t points
+               taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,: ), 'T', -1. )
+      ENDIF
+      !
+   END SUBROUTINE sbc_stress
+
+
    SUBROUTINE sbc_wave( kt )
       !!---------------------------------------------------------------------
@@ -210 +294 @@
          CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing
          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
-      ENDIF
-
-      IF( ln_tauoc .AND. .NOT. cpl_wstrf ) THEN    !==  Wave induced stress  ==!
+         ! check that the drag coefficient contains proper information even if
+         ! the masks do not match - the momentum stress is not masked!
+         WHERE( cdn_wave < 0.0 ) cdn_wave = 1.5e-3
+         WHERE( cdn_wave > 1.0 ) cdn_wave = 1.5e-3
+      ENDIF
+
+      IF( ln_tauoc .AND. .NOT. cpl_tauoc ) THEN    !==  Wave induced stress  ==!
          CALL fld_read( kt, nn_fsbc, sf_tauoc )          ! read wave norm stress from external forcing
          tauoc_wave(:,:) = sf_tauoc(1)%fnow(:,:,1)
-      ENDIF
-
-      IF( ln_sdw )  THEN                           !==  Computation of the 3d Stokes Drift  ==! 
+         WHERE( tauoc_wave <   0.0 ) tauoc_wave = 1.0
+         WHERE( tauoc_wave > 100.0 ) tauoc_wave = 1.0
+      ENDIF
+
+      IF( ln_tauw .AND. .NOT. cpl_tauw ) THEN      !==  Wave induced stress  ==!
+         CALL fld_read( kt, nn_fsbc, sf_tauw )           ! read ocean stress components from external forcing (T grid)
+         tauw_x(:,:) = sf_tauw(1)%fnow(:,:,1)
+         WHERE( tauw_x < -100.0 ) tauw_x = 0.0
+         WHERE( tauw_x >  100.0 ) tauw_x = 0.0
+
+         tauw_y(:,:) = sf_tauw(2)%fnow(:,:,1)
+         WHERE( tauw_y < -100.0 ) tauw_y = 0.0
+         WHERE( tauw_y >  100.0 ) tauw_y = 0.0
+      ENDIF
+
+      IF( ln_phioc .AND. .NOT. cpl_phioc ) THEN    !==  Wave to ocean energy  ==!
+         CALL fld_read( kt, nn_fsbc, sf_phioc )          ! read wave to ocean energy from external forcing
+         rn_crban(:,:) = 29.0 * sf_phioc(1)%fnow(:,:,1)     ! ! Alfa is phioc*sqrt(rau0/zrhoa)  : rau0=water density, zhroa= air density
+         WHERE( rn_crban > 1.e8   ) rn_crban = 0.0    !remove first mask mistmatch points, then cap values in case of low friction velocity
+         WHERE( rn_crban < 0.0    ) rn_crban = 0.0
+         WHERE( rn_crban > 1000.0 ) rn_crban = 1000.0
+      ENDIF
+
+      IF( ln_sdw .OR. ln_rough )  THEN             !==  Computation of the 3d Stokes Drift  ==! 
          !
          IF( jpfld > 0 ) THEN                            ! Read from file only if the field is not coupled
             CALL fld_read( kt, nn_fsbc, sf_sd )          ! read wave parameters from external forcing
-            IF( jp_hsw > 0 )   hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
-            IF( jp_wmp > 0 )   wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
-            IF( jp_usd > 0 )   ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
-            IF( jp_vsd > 0 )   vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
-         ENDIF
-         !
+            IF( jp_hsw > 0 ) THEN
+               hsw  (:,:) = sf_sd(jp_hsw)%fnow(:,:,1)   ! significant wave height
+               WHERE( hsw > 100.0 ) hsw = 0.0
+               WHERE( hsw <   0.0 ) hsw = 0.0
+            ENDIF
+            IF( jp_wmp > 0 ) THEN
+               wmp  (:,:) = sf_sd(jp_wmp)%fnow(:,:,1)   ! wave mean period
+               WHERE( wmp > 100.0 ) wmp = 0.0
+               WHERE( wmp <   0.0 ) wmp = 0.0
+            ENDIF
+            IF( jp_wfr > 0 ) THEN
+               wfreq(:,:) = sf_sd(jp_wfr)%fnow(:,:,1)   ! Peak wave frequency 
+               WHERE( wfreq <    0.0 ) wfreq = 0.0 
+               WHERE( wfreq >  100.0 ) wfreq = 0.0
+            ENDIF
+            IF( jp_usd > 0 ) THEN
+               ut0sd(:,:) = sf_sd(jp_usd)%fnow(:,:,1)   ! 2D zonal Stokes Drift at T point
+               WHERE( ut0sd < -100.0 ) ut0sd = 1.0
+               WHERE( ut0sd >  100.0 ) ut0sd = 1.0
+            ENDIF
+            IF( jp_vsd > 0 ) THEN
+               vt0sd(:,:) = sf_sd(jp_vsd)%fnow(:,:,1)   ! 2D meridional Stokes Drift at T point
+               WHERE( vt0sd < -100.0 ) vt0sd = 1.0
+               WHERE( vt0sd >  100.0 ) vt0sd = 1.0
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      IF( ln_sdw ) THEN
          ! Read also wave number if needed, so that it is available in coupling routines
          IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
@@ -235 +367 @@
          !                                         !==  Computation of the 3d Stokes Drift  ==! 
          !
-         IF( jpfld == 4 )   CALL sbc_stokes()            ! Calculate only if required fields are read
-         !                                               ! In coupled wave model-NEMO case the call is done after coupling
+         IF( ((nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) .AND. &
+                          jp_hsw>0 .AND. jp_wmp>0 .AND. jp_usd>0 .AND. jp_vsd>0) .OR. &
+             (nn_sdrift==jp_peakph .AND. jp_wfr>0 .AND. jp_usd>0 .AND. jp_vsd>0) ) &
+            CALL sbc_stokes()            ! Calculate only if required fields are read
+         !                               ! In coupled wave model-NEMO case the call is done after coupling
          !
       ENDIF
@@ -261 +396 @@
       !!
       CHARACTER(len=100)     ::  cn_dir                          ! Root directory for location of drag coefficient files
-      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i     ! array of namelist informations on the fields to read
-      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd,  &
-                             &   sn_hsw, sn_wmp, sn_wnum, sn_tauoc      ! informations about the fields to be read
-      !
-      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i, slf_j     ! array of namelist informations on the fields to read
+      TYPE(FLD_N)            ::  sn_cdg, sn_usd, sn_vsd, sn_phioc, &
+                             &   sn_hsw, sn_wmp, sn_wfr, sn_wnum , &
+                             &   sn_tauoc, sn_tauwx, sn_tauwy      ! information about the fields to be read
+      !
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wfr, sn_wnum, sn_phioc,          &
+                             sn_tauoc, sn_tauwx, sn_tauwy,                                                       &
+                             ln_cdgw, ln_sdw, ln_stcor, ln_phioc, ln_tauoc, ln_tauw, ln_zdfqiao, ln_rough,       &
+                             nn_sdrift, nn_wmix
       !!---------------------------------------------------------------------
       !
@@ -277 +416 @@
       IF(lwm) WRITE ( numond, namsbc_wave )
       !
+      IF(lwp) THEN               !* Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_wave_init: wave physics'
+         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) '   Namelist namsbc_wave : set wave physics parameters'
+         WRITE(numout,*) '      Stokes drift corr. to vert. velocity    ln_sdw      = ', ln_sdw 
+         WRITE(numout,*) '        vertical parametrization              nn_sdrift   = ', nn_sdrift 
+         WRITE(numout,*) '      Stokes coriolis term                    ln_stcor    = ', ln_stcor 
+         WRITE(numout,*) '      wave modified ocean stress              ln_tauoc    = ', ln_tauoc 
+         WRITE(numout,*) '      wave modified ocean stress components   ln_tauw     = ', ln_tauw 
+         WRITE(numout,*) '      wave to ocean energy                    ln_phioc    = ', ln_phioc
+         WRITE(numout,*) '        vertical mixing parametrization       nn_wmix     = ', nn_wmix 
+         WRITE(numout,*) '      neutral drag coefficient                ln_cdgw     = ', ln_cdgw
+         WRITE(numout,*) '      wave roughness length modification      ln_rough    = ', ln_rough 
+         WRITE(numout,*) '      Qiao vertical mixing formulation        ln_zdfqiao  = ', ln_zdfqiao
+      ENDIF
+
+      IF ( ln_wave ) THEN
+         ! Activated wave physics but no wave physics components activated 
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc &
+                                                                    .OR. ln_rough .OR. ln_zdfqiao) )   THEN  
+             CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauoc=F, ln_tauw=F, ln_stcor=F ', &
+                                                      'ln_phioc=F, ln_rough=F, ln_zdfqiao=F' ) 
+         ELSE
+         IF (ln_stcor .AND. .NOT. ln_sdw) & 
+             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
+         IF ( ln_cdgw .AND. .NOT.(nn_drag==jp_ukmo .OR. nn_drag==jp_std .OR. nn_drag==jp_const .OR. nn_drag==jp_mcore) ) & 
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation must be 0, 1, 2, or 3')
+         IF ( ln_cdgw .AND. ln_blk_core .AND. nn_drag==0 ) &
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in core forcing must be 1, 2, or 3')
+         IF ( ln_cdgw .AND. ln_flx .AND. nn_drag==3 ) &
+             CALL ctl_stop( 'The chosen nn_drag for momentum calculation in direct forcing must be 0, 1, or 2')
+         IF( ln_phioc .AND. .NOT.(nn_wmix==jp_craigbanner .OR. nn_wmix==jp_janssen) ) & 
+            CALL ctl_stop( 'The chosen nn_wmix for wave vertical mixing must be 0, or 1' )
+         IF( ln_sdw .AND. .NOT.(nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips .OR. nn_sdrift==jp_peakph) ) & 
+            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
+         IF( ln_zdfqiao .AND. .NOT.ln_sdw ) & 
+            CALL ctl_stop( 'Qiao vertical mixing can not be used without Stokes drift (ln_sdw)' )
+         IF( ln_tauoc .AND. ln_tauw ) & 
+            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
+                                     '(ln_tauoc=.true. and ln_tauw=.true.)' )
+         IF( ln_tauoc ) &
+             CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauoc=.true.)' )
+         IF( ln_tauw ) &
+             CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', &
+                                  'This will override any other specification of the ocean stress' ) 
+         ENDIF
+      ELSE
+         IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauoc .OR. ln_tauw .OR. ln_stcor .OR. ln_phioc .OR. ln_rough .OR. ln_zdfqiao ) &  
+            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
+            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
+            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 & 
+            &                  'or Stokes-Coriolis term (ln_stcor=T)',                         &
+            &                  'or ocean stress modification due to waves (ln_tauoc=T) ',      & 
+            &                  'or ocean stress components from waves (ln_tauw=T) ',          & 
+            &                  'or wave to ocean energy modification (ln_phioc=T) ',           & 
+            &                  'or wave surface roughness (ln_rough=T) ',                      & 
+            &                  'or Qiao vertical mixing formulation (ln_zdfqiao=T) ' )
+      ENDIF 
+      !
       IF( ln_cdgw ) THEN
          IF( .NOT. cpl_wdrag ) THEN
             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_cd structure' )
             !
                                    ALLOCATE( sf_cd(1)%fnow(jpi,jpj,1)   )
             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
+            CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( cdn_wave(jpi,jpj) )
@@ -290 +489 @@
 
       IF( ln_tauoc ) THEN
-         IF( .NOT. cpl_wstrf ) THEN
+         IF( .NOT. cpl_tauoc ) THEN
             ALLOCATE( sf_tauoc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_tauoc
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauoc structure' )
             !
                                     ALLOCATE( sf_tauoc(1)%fnow(jpi,jpj,1)   )
             IF( sn_tauoc%ln_tint )  ALLOCATE( sf_tauoc(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'Wave module', 'namsbc_wave' )
+            CALL fld_fill( sf_tauoc, (/ sn_tauoc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( tauoc_wave(jpi,jpj) )
       ENDIF
 
-      IF( ln_sdw ) THEN   ! Find out how many fields have to be read from file if not coupled
-         jpfld=0
-         jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0
-         IF( .NOT. cpl_sdrftx ) THEN
+      IF( ln_tauw ) THEN
+         IF( .NOT. cpl_tauw ) THEN
+            ALLOCATE( sf_tauw(2), STAT=ierror )           !* allocate and fill sf_wave with sn_tauwx/y
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_tauw structure' )
+            !
+            ALLOCATE( slf_j(2) )
+            slf_j(1) = sn_tauwx
+            slf_j(2) = sn_tauwy
+                                    ALLOCATE( sf_tauw(1)%fnow(jpi,jpj,1)   )
+                                    ALLOCATE( sf_tauw(2)%fnow(jpi,jpj,1)   )
+            IF( slf_j(1)%ln_tint )  ALLOCATE( sf_tauw(1)%fdta(jpi,jpj,1,2) )
+            IF( slf_j(2)%ln_tint )  ALLOCATE( sf_tauw(2)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill( sf_tauw, (/ slf_j /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( tauw_x(jpi,jpj) )
+         ALLOCATE( tauw_y(jpi,jpj) )
+      ENDIF
+
+      IF( ln_phioc ) THEN
+         IF( .NOT. cpl_phioc ) THEN
+            ALLOCATE( sf_phioc(1), STAT=ierror )           !* allocate and fill sf_wave with sn_phioc
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_phioc structure' )
+            !
+                                    ALLOCATE( sf_phioc(1)%fnow(jpi,jpj,1)   )
+            IF( sn_phioc%ln_tint )  ALLOCATE( sf_phioc(1)%fdta(jpi,jpj,1,2) )
+            CALL fld_fill( sf_phioc, (/ sn_phioc /), cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+         ENDIF
+         ALLOCATE( rn_crban(jpi,jpj) )
+      ENDIF
+
+      ! Find out how many fields have to be read from file if not coupled
+      jpfld=0
+      jp_usd=0   ;   jp_vsd=0   ;   jp_hsw=0   ;   jp_wmp=0   ;   jp_wfr=0
+      IF( ln_sdw ) THEN
+         IF( .NOT. cpl_sdrft ) THEN
             jpfld  = jpfld + 1
             jp_usd = jpfld
-         ENDIF
-         IF( .NOT. cpl_sdrfty ) THEN
             jpfld  = jpfld + 1
             jp_vsd = jpfld
          ENDIF
-         IF( .NOT. cpl_hsig ) THEN
+         IF( .NOT. cpl_hsig .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
             jpfld  = jpfld + 1
             jp_hsw = jpfld
          ENDIF
-         IF( .NOT. cpl_wper ) THEN
+         IF( .NOT. cpl_wper .AND. (nn_sdrift==jp_breivik .OR. nn_sdrift==jp_phillips) ) THEN
             jpfld  = jpfld + 1
             jp_wmp = jpfld
          ENDIF
-
-         ! Read from file only the non-coupled fields 
-         IF( jpfld > 0 ) THEN
-            ALLOCATE( slf_i(jpfld) )
-            IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
-            IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
-            IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
-            IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
-            ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_wave structure' )
-            !
-            DO ifpr= 1, jpfld
-               ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
-               IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
-            END DO
-            !
-            CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'Wave module ', 'namsbc_wave' )
-         ENDIF
+         IF( .NOT. cpl_wfreq .AND. nn_sdrift==jp_peakph ) THEN
+            jpfld  = jpfld + 1
+            jp_wfr = jpfld
+         ENDIF
+      ENDIF
+
+      IF( ln_rough .AND. .NOT. cpl_hsig .AND. jp_hsw==0 ) THEN
+         jpfld  = jpfld + 1
+         jp_hsw = jpfld
+      ENDIF
+
+      ! Read from file only the non-coupled fields 
+      IF( jpfld > 0 ) THEN
+         ALLOCATE( slf_i(jpfld) )
+         IF( jp_usd > 0 )   slf_i(jp_usd) = sn_usd
+         IF( jp_vsd > 0 )   slf_i(jp_vsd) = sn_vsd
+         IF( jp_hsw > 0 )   slf_i(jp_hsw) = sn_hsw
+         IF( jp_wmp > 0 )   slf_i(jp_wmp) = sn_wmp
+         IF( jp_wfr > 0 )   slf_i(jp_wfr) = sn_wfr
+         ALLOCATE( sf_sd(jpfld), STAT=ierror )   !* allocate and fill sf_sd with stokes drift
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable to allocate sf_sd structure' )
+         !
+         DO ifpr= 1, jpfld
+            ALLOCATE( sf_sd(ifpr)%fnow(jpi,jpj,1) )
+            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf_sd(ifpr)%fdta(jpi,jpj,1,2) )
+         END DO
+         !
+         CALL fld_fill( sf_sd, slf_i, cn_dir, 'sbc_wave_init', 'read wave input', 'namsbc_wave' )
+      ENDIF
+
+      IF( ln_sdw ) THEN
          ALLOCATE( usd  (jpi,jpj,jpk), vsd  (jpi,jpj,jpk), wsd(jpi,jpj,jpk) )
-         ALLOCATE( hsw  (jpi,jpj)    , wmp  (jpi,jpj)     )
+         ALLOCATE( wmp  (jpi,jpj)  )
+         ALLOCATE( wfreq (jpi,jpj) )
          ALLOCATE( ut0sd(jpi,jpj)    , vt0sd(jpi,jpj)     )
          ALLOCATE( div_sd(jpi,jpj) )
@@ -349 +591 @@
          IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
             ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
-            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wave structure' )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_init: unable toallocate sf_wn structure' )
                                    ALLOCATE( sf_wn(1)%fnow(jpi,jpj,1)   )
             IF( sn_wnum%ln_tint )  ALLOCATE( sf_wn(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'Wave module', 'namsbc_wave' )
+            CALL fld_fill( sf_wn, (/ sn_wnum /), cn_dir, 'sbc_wave', 'read wave input', 'namsbc_wave' )
          ENDIF
          ALLOCATE( wnum(jpi,jpj) )
+      ENDIF
+
+      IF( ln_sdw .OR. ln_rough ) THEN
+         ALLOCATE( hsw  (jpi,jpj) )
       ENDIF
       !

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

@@ -10 +10 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
+   USE sbc_oce, ONLY : ln_zdfqiao
 
    IMPLICIT NONE
@@ -35 +36 @@
    INTEGER , PUBLIC ::   nn_npc      !: non penetrative convective scheme call  frequency
    INTEGER , PUBLIC ::   nn_npcp     !: non penetrative convective scheme print frequency
-   LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag
 
 

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

@@ -24 +24 @@
    USE phycst         ! physical constants
    USE zdfmxl         ! mixed layer
-   USE sbcwave ,  ONLY: hsw   ! significant wave height  
+   USE sbcwave
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -48 +48 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustarb2 !: Squared bottom  velocity scale at T-points
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke1 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_tke3 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rsbc_psi1 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rtrans
 
    !                              !! ** Namelist  namzdf_gls  **
@@ -61 +66 @@
    REAL(wp) ::   rn_emin           ! minimum value of TKE (m2/s2)
    REAL(wp) ::   rn_charn          ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value)
-   REAL(wp) ::   rn_crban          ! Craig and Banner constant for surface breaking waves mixing
+   REAL(wp) ::   rn_crban_default  ! Craig and Banner constant for surface breaking waves mixing
    REAL(wp) ::   rn_hsro           ! Minimum surface roughness
-   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 
+   REAL(wp) ::   rn_frac_hs        ! Fraction of wave height as surface roughness (if nn_z0_met = 1, 3) 
 
    REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters
@@ -92 +97 @@
    REAL(wp) ::   rm7           =  0.0_wp
    REAL(wp) ::   rm8           =  0.318_wp
-   REAL(wp) ::   rtrans        =  0.1_wp
+   REAL(wp) ::   rtrans_default =  0.1_wp
    REAL(wp) ::   rc02, rc02r, rc03, rc04                          ! coefficients deduced from above parameters
-   REAL(wp) ::   rsbc_tke1, rsbc_tke2, rfact_tke                  !     -           -           -        -
-   REAL(wp) ::   rsbc_psi1, rsbc_psi2, rfact_psi                  !     -           -           -        -
+   REAL(wp) ::   rsbc_tke1_default, rsbc_tke2, rfact_tke          !     -           -           -        - 
+   REAL(wp) ::   rsbc_psi2, rsbc_psi3, rfact_psi                  !     -           -           -        -
    REAL(wp) ::   rsbc_zs1, rsbc_zs2                               !     -           -           -        -
    REAL(wp) ::   rc0, rc2, rc3, rf6, rcff, rc_diff                !     -           -           -        -
@@ -117 +122 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
-         &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
+      ALLOCATE( mxln(jpi,jpj,jpk) , zwall(jpi,jpj,jpk) ,     & 
+         &      ustars2(jpi,jpj)  , ustarb2(jpi,jpj)   ,     & 
+         &      rsbc_tke1(jpi,jpj), rsbc_tke3(jpi,jpj) ,     & 
+         &      rsbc_psi1(jpi,jpj), rtrans(jpi,jpj), STAT= zdf_gls_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_gls_alloc )
@@ -163 +170 @@
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
+      ! variable initialization 
+      IF( ln_phioc ) THEN 
+         IF( nn_wmix==jp_janssen ) THEN 
+            rsbc_tke1(:,:) = (-rsc_tke*rn_crban(:,:)/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp)  ! k_eps = 53.Dirichlet + Wave breaking  
+            rsbc_tke3(:,:) = rdt * rn_crban(:,:)                                           ! Neumann + Wave breaking 
+            rsbc_psi1(:,:) = rc0**rpp * rsbc_tke1(:,:)**rmm * rl_sf**rnn                   ! Dirichlet + Wave breaking 
+         ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+            rsbc_tke1(:,:) = -3._wp/2._wp*rn_crban(:,:)*ra_sf*rl_sf 
+            rsbc_tke3(:,:) = rdt * rn_crban(:,:) / rl_sf 
+            rtrans(:,:) = 0.2_wp / MAX( rsmall, rsbc_tke1(:,:)**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
+         ENDIF 
+      ENDIF 
+
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
          avt (:,:,:) = avt_k (:,:,:)
@@ -200 +220 @@
          zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)  
-         zhsro(:,:) = hsw(:,:)
+         zhsro = MAX(rn_frac_hs * hsw, rn_hsro)
       END SELECT
 
@@ -314 +334 @@
       ! ---------------------------------
       !
+      IF( ln_phioc ) THEN 
+         SELECT CASE ( nn_bc_surf ) 
+         ! 
+         CASE ( 0 )             ! Dirichlet case 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! First level 
+               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+ 
+               ! One level below 
+               en(:,:,2) = MAX( rsbc_tke1(:,:) * ustars2(:,:) * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 
+               en(:,:,1) = MAX(en(:,:,1), rn_emin)  
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               !  
+               ! One level below 
+               en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:) * ((zhsro(:,:)+fsdepw(:,:,2)) & 
+                   &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 
+               en(:,:,2) = MAX(en(:,:,2), rn_emin ) 
+               z_elem_a(:,:,2) = 0._wp  
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+               ! 
+            ENDIF 
+         CASE ( 1 )             ! Neumann boundary condition on d(e)/dz 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! Dirichlet conditions at k=1 
+               en(:,:,1) = MAX( rsbc_tke1(:,:) * ustars2(:,:), rn_emin ) 
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! at k=2, set de/dz=Fw 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               zflxs(:,:) = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * ((zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 
+               en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               ! Dirichlet conditions at k=1 
+               en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1(:,:))**(2._wp/3._wp) 
+               en(:,:,1)       = MAX(en(:,:,1), rn_emin)       
+               z_elem_a(:,:,1) = en(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! at k=2, set de/dz=Fw 
+               !cbr 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:)) )) 
+               zflxs(:,:)      = rsbc_tke3(:,:) * ustars2(:,:)**1.5_wp * zkar(:,:) & 
+                    &                           * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 
+
+               en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 
+               ! 
+            ENDIF 
+         END SELECT 
+      ELSE
       SELECT CASE ( nn_bc_surf )
       !
       CASE ( 0 )             ! Dirichlet case
       ! First level
-      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
       en(:,:,1) = MAX(en(:,:,1), rn_emin) 
       z_elem_a(:,:,1) = en(:,:,1)
@@ -325 +410 @@
       ! 
       ! One level below
-      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) &
+      en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default * ((zhsro(:,:)+fsdepw(:,:,2)) &
           &            / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
       en(:,:,2) = MAX(en(:,:,2), rn_emin )
@@ -334 +419 @@
       !
       CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
-      !
       ! Dirichlet conditions at k=1
-      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp)
+      en(:,:,1)       = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1_default)**(2._wp/3._wp)
       en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
       z_elem_a(:,:,1) = en(:,:,1)
@@ -346 +430 @@
       z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
       z_elem_a(:,:,2) = 0._wp
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) ))
+      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:)) ))
       zflxs(:,:)      = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) &
            &                      * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf)
@@ -354 +438 @@
       !
       END SELECT
+      ENDIF ! ln_phioc
 
       ! Bottom boundary condition on tke
@@ -539 +624 @@
       !
       CASE ( 0 )             ! Dirichlet boundary conditions
-      !
-      ! Surface value
-      zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
-      !
-      ! One level below
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) )))
-      zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)
-      psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,2) = 0._wp
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
-      ! 
-      !
+         IF( ln_phioc ) THEN   ! Wave induced mixing case 
+                               ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 
+                               ! balance between the production and the 
+                               ! dissipation terms including the wave effect 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! Surface value 
+               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1)       
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! One level below 
+               zex1 = (rmm*ra_sf+rnn) 
+               zex2 = (rmm*ra_sf) 
+               zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 
+               psi (:,:,2) = rsbc_psi1(:,:) * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+               !  
+               ! 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               ! Surface value 
+               zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! One level below 
+               zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdepw(:,:,2)/zhsro(:,:) ))) 
+               zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
+               psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,2) = 0._wp 
+               z_elem_c(:,:,2) = 0._wp 
+               z_elem_b(:,:,2) = 1._wp 
+               !  
+               ! 
+            ENDIF 
+         ELSE                  ! Wave induced mixing case with default values 
+            ! Surface value 
+            zdep(:,:)       = zhsro(:,:) * rl_sf ! Cosmetic 
+            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+            z_elem_a(:,:,1) = psi(:,:,1) 
+            z_elem_c(:,:,1) = 0._wp 
+            z_elem_b(:,:,1) = 1._wp 
+            ! 
+            ! One level below 
+            zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdepw(:,:,2)/zhsro(:,:) ))) 
+            zdep(:,:)       = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 
+            psi (:,:,2)     = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+            z_elem_a(:,:,2) = 0._wp 
+            z_elem_c(:,:,2) = 0._wp 
+            z_elem_b(:,:,2) = 1._wp 
+            !  
+            ! 
+         ENDIF
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
-      !
-      ! Surface value: Dirichlet
-      zdep(:,:)       = zhsro(:,:) * rl_sf
-      psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-      z_elem_a(:,:,1) = psi(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
-      !
-      ! Neumann condition at k=2
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
-      !
-      ! Set psi vertical flux at the surface:
-      zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
-      zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
-      zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
-      zdep(:,:) =  rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * &
-             & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)
-      zflxs(:,:) = zdep(:,:) * zflxs(:,:)
-      psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)
-
-      !   
-      !
+         IF( ln_phioc ) THEN  ! Wave induced mixing case with forced/coupled fields 
+            IF( nn_wmix==jp_janssen ) THEN 
+               ! 
+               zdep(:,:) = rl_sf * zhsro(:,:) 
+               psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! Neumann condition at k=2 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               ! 
+               ! Set psi vertical flux at the surface: 
+               zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 
+               zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 
+                  &                   * en(:,:,1)**rmm * zdep 
+               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+               ! 
+            ELSE IF( nn_wmix==jp_craigbanner ) THEN 
+               ! Surface value: Dirichlet 
+               zdep(:,:)       = zhsro(:,:) * rl_sf 
+               psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+               z_elem_a(:,:,1) = psi(:,:,1) 
+               z_elem_c(:,:,1) = 0._wp 
+               z_elem_b(:,:,1) = 1._wp 
+               ! 
+               ! Neumann condition at k=2 
+               z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+               z_elem_a(:,:,2) = 0._wp 
+               ! 
+               ! Set psi vertical flux at the surface: 
+               zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans(:,:)*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
+               zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
+               zflxs(:,:) = (rnn + rsbc_tke1(:,:) * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
+                          (1._wp + rsbc_tke1(:,:) * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
+               zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
+                      & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
+               zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
+               psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+               !    
+            ENDIF 
+         ELSE                 ! Wave induced mixing case with default values 
+            ! Surface value: Dirichlet 
+            zdep(:,:)       = zhsro(:,:) * rl_sf 
+            psi (:,:,1)     = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 
+            z_elem_a(:,:,1) = psi(:,:,1) 
+            z_elem_c(:,:,1) = 0._wp 
+            z_elem_b(:,:,1) = 1._wp 
+            ! 
+            ! Neumann condition at k=2 
+            z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 
+            z_elem_a(:,:,2) = 0._wp 
+            ! 
+            ! Set psi vertical flux at the surface: 
+            zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans_default*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 
+            zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 
+            zflxs(:,:) = (rnn + rsbc_tke1_default * (rnn + rmm*ra_sf) * zdep(:,:)) * & 
+                       (1._wp + rsbc_tke1_default * zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 
+            zdep(:,:) =  rsbc_psi1(:,:) * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 
+                   & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 
+            zflxs(:,:) = zdep(:,:) * zflxs(:,:) 
+            psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 
+            !    
+            ! 
+         ENDIF
       END SELECT
 
@@ -870 +1039 @@
       NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, &
          &            rn_clim_galp, ln_sigpsi, rn_hsro,      &
-         &            rn_crban, rn_charn, rn_frac_hs,        &
+         &            rn_crban_default, rn_charn, rn_frac_hs,&
          &            nn_bc_surf, nn_bc_bot, nn_z0_met,      &
          &            nn_stab_func, nn_clos
@@ -898 +1067 @@
          WRITE(numout,*) '      TKE Bottom boundary condition                 nn_bc_bot      = ', nn_bc_bot
          WRITE(numout,*) '      Modify psi Schmidt number (wb case)           ln_sigpsi      = ', ln_sigpsi
-         WRITE(numout,*) '      Craig and Banner coefficient                  rn_crban       = ', rn_crban
+         WRITE(numout,*) '      Craig and Banner coefficient (default)        rn_crban       = ', rn_crban_default
          WRITE(numout,*) '      Charnock coefficient                          rn_charn       = ', rn_charn
          WRITE(numout,*) '      Surface roughness formula                     nn_z0_met      = ', nn_z0_met
@@ -915 +1084 @@
       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' )  
       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' )  
-      IF( nn_z0_met == 3 .AND. .NOT.ln_wave ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_wave=T' )  
+      IF( nn_z0_met == 3 .AND. .NOT.ln_rough ) CALL ctl_stop( 'zdf_gls_init: nn_z0_met=3 requires ln_rough=T' )  
+      IF( nn_z0_met .NE. 3 .AND. ln_rough ) THEN 
+         CALL ctl_warn('W A R N I N G:  ln_rough=.TRUE. but nn_z0_met is not 3 - resetting nn_z0_met to 3')  
+         nn_z0_met = 3  
+      ENDIF
       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' )  
       IF( nn_clos       < 0 .OR. nn_clos       > 3 ) CALL ctl_stop( 'zdf_gls_init: bad flag: nn_clos is 0, 1, 2 or 3' )
@@ -1089 +1262 @@
                &                              - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) )
 
-      IF ( rn_crban==0._wp ) THEN
+      IF( .NOT. ln_phioc .AND. rn_crban_default==0._wp ) THEN
          rl_sf = vkarmn
       ELSE
@@ -1128 +1301 @@
       rc03  = rc02 * rc0
       rc04  = rc03 * rc0
-      rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf                      ! Dirichlet + Wave breaking
-      rsbc_tke2 = rdt * rn_crban / rl_sf                                 ! Neumann + Wave breaking 
-      zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp )
-      rtrans = 0.2_wp / zcr                                              ! Ad. inverse transition length between log and wave layer 
+      rsbc_tke1_default = -3._wp/2._wp*rn_crban_default*ra_sf*rl_sf 
+      rsbc_tke2         = rdt * rn_crban_default / rl_sf 
+      zcr               = MAX( rsmall, rsbc_tke1_default**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 
+      rtrans_default    = 0.2_wp / zcr
       rsbc_zs1  = rn_charn/grav                                          ! Charnock formula for surface roughness
       rsbc_zs2  = rn_frac_hs / 0.85_wp / grav * 665._wp                  ! Rascle formula for surface roughness 

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

@@ -54 +54 @@
       !!
       NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  &  
-         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp,       &  
-         &        ln_zdfqiao
+         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp
       !!----------------------------------------------------------------------
 
@@ -84 +83 @@
          WRITE(numout,*) '      npc call  frequency                 nn_npc    = ', nn_npc
          WRITE(numout,*) '      npc print frequency                 nn_npcp   = ', nn_npcp
-         WRITE(numout,*) '      Qiao formulation flag               ln_zdfqiao=', ln_zdfqiao
       ENDIF
 

branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -106 +106 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )  THEN
-         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
-                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-      ENDIF
-                         CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+                         CALL sbc     ( kstp )        ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>