Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7716)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7717)
@@ -299,5 +299,6 @@
    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 
 /
 !-----------------------------------------------------------------------
@@ -1296,5 +1297,8 @@
    sn_tauoc    =  'sdw_wave' ,        1          , 'wave_stress',     .true.   , .false. , 'daily'   ,  ''      , ''       , '' 
 ! 
-   cn_dir  = './'  !  root directory for the location of drag coefficient files
+   cn_dir      = './'  !  root directory for the location of drag coefficient files
+   nn_drag     = 0     !  formula to calculate momentum from the wind components
+                              =0 UKMO SHELF formulation
+                              =1 standard formulation
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7716)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7717)
@@ -956,4 +956,5 @@
       !!----------------------------------------------------------------------
       USE zdf_oce,  ONLY : ln_zdfqiao
+      USE sbcflx ,  ONLY : ln_shelf_flx
 
       INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
@@ -1039,7 +1040,22 @@
       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, jpj
+               DO ji = 1, jpi
+                  ! here utau and vtau should contain the wind components as read from the forcing files
+                  zcoef = SQRT(utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj))
+                  frcv(jpr_otx1)%z3(ji,jj,1) = zrhoa * frcv(jpr_wdrag)%z3(ji,jj,1) * utau(ji,jj) * zcoef
+                  frcv(jpr_oty1)%z3(ji,jj,1) = zrhoa * frcv(jpr_wdrag)%z3(ji,jj,1) * vtau(ji,jj) * zcoef
+                  utau(ji,jj) = frcv(jpr_otx1)%z3(ji,jj,1)
+                  vtau(ji,jj) = frcv(jpr_oty1)%z3(ji,jj,1)
+               END DO
+            END DO
+            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
@@ -1061,4 +1077,7 @@
             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
@@ -1075,4 +1094,5 @@
       !                                                      ! ========================= !
       !                                                      !      10 m wind speed      !   (wndm)
+      !                                                      !   include wave drag coef  !   (wndm)
       !                                                      ! ========================= !
       !
@@ -1085,5 +1105,9 @@
 !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
@@ -1095,4 +1119,5 @@
       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
             IF( srcv(jpr_otx1)%laction ) THEN
@@ -1100,8 +1125,7 @@
                vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
             ENDIF
-            IF( srcv(jpr_taum)%laction .OR. (.NOT.srcv(jpr_taum)%laction .AND. srcv(jpr_otx1)%laction .AND. llnewtx ))   &
+            IF( srcv(jpr_taum)%laction )   &
                taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
-            IF( srcv(jpr_w10m)%laction .OR. (.NOT.srcv(jpr_w10m)%laction .AND. llnewtau .AND.                            &
-               (srcv(jpr_taum)%laction .OR. (.NOT.srcv(jpr_taum)%laction .AND. srcv(jpr_otx1)%laction .AND. llnewtx )))) &
+            IF( srcv(jpr_w10m)%laction )   &
                wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
          ELSE IF( ll_purecpl ) THEN
@@ -1172,9 +1196,4 @@
       IF( srcv(jpr_wstrf)%laction .AND. ln_tauoc ) tauoc_wave(:,:) = frcv(jpr_wstrf)%z3(:,:,1) 
       
-      !                                                      ! ========================= !  
-      !                                                      !   Wave drag coefficient   ! 
-      !                                                      ! ========================= !  
-      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw ) cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
-
       !  Fields received by SAS when OASIS coupling
       !  (arrays no more filled at sbcssm stage)
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 7716)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 7717)
@@ -22,4 +22,5 @@
    USE lib_mpp         ! distribued memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE sbcwave         ! wave physics
 
    IMPLICIT NONE
@@ -35,4 +36,6 @@
    INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
+   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag
+   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx) 
 
    !! * Substitutions
@@ -81,7 +84,9 @@
       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
       !!
       CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files
+      NAMELIST/namsbc_flx/ ln_shelf_flx                           ! Put here to allow merging with another UKMO branch
       TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures
       TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  ! informations about the fields to be read
@@ -140,4 +145,26 @@
             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 .AND. 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) = sf(jp_utau)%fnow(ji,jj,1)
+                     vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
+                  END DO
+               END DO
+            ELSE
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     totwind = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0)
+                     utau(ji,jj) = zrhoa * cdn_wave(ji,jj) * sf(jp_utau)%fnow(ji,jj,1) * totwind
+                     vtau(ji,jj) = zrhoa * cdn_wave(ji,jj) * sf(jp_vtau)%fnow(ji,jj,1) * totwind
+                  END DO
+               END DO
+            ENDIF
+         ENDIF
          !                                                        ! add to qns the heat due to e-p
          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
@@ -153,5 +180,11 @@
                zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
                taum(ji,jj) = zmod
+               IF( ln_shelf_flx ) THEN
+                  ztx = sf(jp_utau)%fnow(ji-1,jj  ) + sf(jp_utau)%fnow(ji,jj)
+                  zty = sf(jp_vtau)%fnow(ji  ,jj-1) + sf(jp_vtau)%fnow(ji,jj)
+                  wndm(ji,jj) = 0.5 * SQRT( ztx * ztx + zty * zty )
+               ELSE
                wndm(ji,jj) = SQRT( zmod * zcoef )
+               ENDIF
             END DO
          END DO
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7716)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 7717)
@@ -137,5 +137,5 @@
          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,*) '                 neutral drag coefficient                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
@@ -225,7 +225,4 @@
          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')
Index: /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7716)
+++ /branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7717)
@@ -63,4 +63,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd        !: surface Stokes drift velocities at t-point
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd  , vsd  , wsd   !: Stokes drift velocities at u-, v- & w-points, resp.
+
+   INTEGER, PUBLIC            ::   nn_drag       ! type of formula to calculate wind stress from wind components
+   INTEGER, PUBLIC, PARAMETER ::   jp_ukmo = 0        ! UKMO SHELF formulation
+   INTEGER, PUBLIC, PARAMETER ::   jp_std  = 1        ! standard formulation: 
 
    !! * Substitutions
@@ -210,4 +214,6 @@
          CALL fld_read( kt, nn_fsbc, sf_cd )             ! read from external forcing
          cdn_wave(:,:) = sf_cd(1)%fnow(:,:,1)
+         ! make sure that the interpolation does not write silly fill values in case it fails in a grid point
+         where( cdn_wave > 1.0 ) cdn_wave = 1.5e-3 
       ENDIF
 
@@ -265,5 +271,5 @@
                              &   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
+      NAMELIST/namsbc_wave/  sn_cdg, cn_dir, sn_usd, sn_vsd, sn_hsw, sn_wmp, sn_wnum, sn_tauoc, nn_drag
       !!---------------------------------------------------------------------
       !
@@ -278,4 +284,10 @@
       !
       IF( ln_cdgw ) THEN
+         ! check the value of nn_drag
+         IF( .NOT. ( nn_drag == jp_ukmo .OR. nn_drag == jp_std )) THEN
+            CALL ctl_stop( 'sbc_wave_init: nn_drag does not specify a valid drag formulation' )
+            RETURN  
+         ENDIF
+           
          IF( .NOT. cpl_wdrag ) THEN
             ALLOCATE( sf_cd(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
