Index: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 12611)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 12654)
@@ -92,5 +92,5 @@
       REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface
       REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector
-      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
+      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at U,V-points
 
       REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3
@@ -167,5 +167,5 @@
          ENDIF
 
-         !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch
+         !! UKMO SHELF flux files contain wind speed not wind stress
          IF( ln_shelf_flx ) THEN
             CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j )
@@ -176,5 +176,12 @@
                   !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing
                   pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1)
-                  !! UKMO SHELF flux files contain wind speed not wind stress
+
+                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
+                  !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
+                  qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1)
+                  !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
+                  emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1)
+
+                  !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch
                   IF( ln_rel_wind ) THEN
                      zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj)
@@ -184,7 +191,19 @@
                      zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
                   ENDIF
-                  wndm(ji,jj) = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj))
-
-                  ! add modification due to drag coefficient read from wave forcing 
+               END DO
+            END DO
+
+            ! Calculate wind speed from the wind components 
+            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
+
+            ! add modification due to drag coefficient read from wave forcing 
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
                   IF( ln_cdgw .AND. nn_drag == jp_std ) THEN 
                      IF( cpl_wdrag ) THEN  
@@ -194,29 +213,18 @@
                         vtau(ji,jj) = zwnd_j(ji,jj) 
                      ELSE 
-                        utau(ji,jj) = zrhoa * cdn_wave(ji,jj) * zwnd_i(ji,jj) * wndm(ji,jj)
-                        vtau(ji,jj) = zrhoa * cdn_wave(ji,jj) * zwnd_j(ji,jj) * wndm(ji,jj)
+                        utau(ji,jj) = zrhoa * 0.25 * ( cdn_wave(ji,jj) + cdn_wave(ji+1,jj) ) * zwnd_i(ji,jj) * &
+                                                     ( wndm(ji,jj)     + wndm(ji+1,jj) )
+                        vtau(ji,jj) = zrhoa * 0.25 * ( cdn_wave(ji,jj) + cdn_wave(ji,jj+1) ) * zwnd_j(ji,jj) * &
+                                                     ( wndm(ji,jj)     + wndm(ji,jj+1) )
                      ENDIF 
                   ELSE IF( nn_drag == jp_const ) THEN 
-                     utau(ji,jj) = zrhoa * zcdrag * zwnd_i(ji,jj) * wndm(ji,jj)
-                     vtau(ji,jj) = zrhoa * zcdrag * zwnd_j(ji,jj) * wndm(ji,jj)
+                     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) )
                   ELSE IF( nn_drag == jp_ukmo ) THEN 
-                     cs = 0.63 + (0.066 * wndm(ji,jj))
-                     utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * wndm(ji,jj)
-                     vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * wndm(ji,jj)
+                     totwindspd = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj))
+                     cs = 0.63 + (0.066 * totwindspd)
+                     utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * totwindspd
+                     vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * totwindspd
                   ENDIF 
-                  taum(ji,jj) = sqrt(utau(ji,jj)*utau(ji,jj) + vtau(ji,jj)*vtau(ji,jj))
-
-                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
-                  !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot)
-                  qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1)
-                  !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P
-                  emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1)
-               END DO
-            END DO
-            ! Move tau components to the right grid
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  utau(ji,jj) = 0.5 * ( utau(ji,jj) + utau(ji+1,jj) )
-                  vtau(ji,jj) = 0.5 * ( vtau(ji,jj) + vtau(ji,jj+1) )
 
                   !!UKMO SHELF effect of atmospheric pressure on SSH
@@ -228,4 +236,6 @@
                END DO
             END DO
+
+            CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j )
          ELSE
             DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
@@ -245,18 +255,18 @@
                END DO
             END DO
+         ENDIF
             !                                                        ! module of wind stress and wind speed at T-point
-            zcoef = 1. / ( zrhoa * zcdrag )
+         zcoef = 1. / ( zrhoa * zcdrag )
 !CDIR NOVERRCHK
-            DO jj = 2, jpjm1
+         DO jj = 2, jpjm1
 !CDIR NOVERRCHK
-               DO ji = fs_2, fs_jpim1   ! vect. opt.
-                  ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-                  zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-                  zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-                  taum(ji,jj) = zmod
-                  wndm(ji,jj) = SQRT( zmod * zcoef )
-               END DO
-            END DO
-         ENDIF
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               ztx = utau(ji-1,jj  ) + utau(ji,jj) 
+               zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
+               zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+               taum(ji,jj) = zmod
+               IF( .NOT. ln_shelf_flx ) wndm(ji,jj) = SQRT( zmod * zcoef )
+            END DO
+         END DO
          !                                                        ! add to qns the heat due to e-p
          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
@@ -284,8 +294,4 @@
          ENDIF
          !
-         IF( ln_shelf_flx ) THEN
-            CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j )
-         ENDIF
-         !
       ENDIF
       !
