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 12598)
+++ /branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 12599)
@@ -135,14 +135,14 @@
          slf_i(jp_emp ) = sn_emp
          !
-            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
-            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press
-   
-            ! define local jpfld depending on shelf_flx logical
-            IF( ln_shelf_flx ) THEN
-               jpfld_local = jpfld
-            ELSE
-               jpfld_local = jpfld-1
-            ENDIF
-            !
+         ! define local jpfld depending on shelf_flx logical
+         IF( ln_shelf_flx ) THEN
+            slf_i(jp_press) = sn_press
+            jpfld_local = jpfld
+         ELSE
+            jpfld_local = jpfld-1
+         ENDIF
+         !
+         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure
+         !
          IF( ierror > 0 ) THEN   
             CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN  
@@ -163,133 +163,110 @@
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency
 
+         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
+         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
+         ENDIF
+
          !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch
          IF( ln_shelf_flx ) THEN
             CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j )
 
-            IF( ln_rel_wind ) THEN
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
+            ! set the ocean fluxes from read fields
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !! 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
+                  IF( ln_rel_wind ) THEN
                      zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj)
                      zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) - rn_wfac * ssv_m(ji,jj)
-                  END DO
-               END DO
-            ELSE
-               zwnd_i(:,:) = sf(jp_utau)%fnow(:,:,1)
-               zwnd_j(:,:) = sf(jp_vtau)%fnow(:,:,1)
-            ENDIF
-         ENDIF
-
-         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle
-         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1)
-         ENDIF
-!CDIR COLLAPSE
-            !!UKMO SHELF effect of atmospheric pressure on SSH
-            ! If using ln_apr_dyn, this is done there so don't repeat here.
-            IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
+                  ELSE
+                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
+                     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 
+                  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 
+                        utau(ji,jj) = zwnd_i(ji,jj) 
+                        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)
+                     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)
+                  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)
+                  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
+                  ! If using ln_apr_dyn, this is done there so don't repeat here.
+                  IF( .NOT. ln_apr_dyn) THEN
                      apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj)
                      apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj)
-                  END DO
-               END DO
-            ENDIF ! ln_shelf_flx
-
-         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
-            DO ji = 1, jpi
-                   IF( ln_shelf_flx ) THEN
-                      !! 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
-                      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
-                   ELSE
-                      utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
-                      vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
-                   ENDIF
-                   qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
-                   IF( ln_foam_flx .OR. ln_shelf_flx ) THEN
-                      !! 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)
-                   ELSE
-                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
-                      emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
-                   ENDIF
-            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
+               END DO
+            END DO
+         ELSE
+            DO jj = 1, jpj                                           ! set the ocean fluxes from read fields
+               DO ji = 1, jpi
+                  utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1)
+                  vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1)
+                  qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1)
+                  IF( ln_foam_flx ) THEN
+                     !! 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)
+                  ELSE
+                     qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1)
+                     emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1)
+                  ENDIF
+               END DO
+            END DO
+            !                                                        ! module of wind stress and wind speed at T-point
+            zcoef = 1. / ( zrhoa * zcdrag )
+!CDIR NOVERRCHK
+            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
          !                                                        ! add to qns the heat due to e-p
          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST
          !
-   
-            !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe
-            IF( ln_foam_flx ) THEN
-               CALL lbc_lnk( utau(:,:), 'U', -1. )
-               CALL lbc_lnk( vtau(:,:), 'V', -1. )
-            ENDIF
+         !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe
+         IF( ln_foam_flx .OR. ln_shelf_flx ) THEN
+            CALL lbc_lnk_multi( utau(:,:), 'U', -1., vtau(:,:), 'V', -1. ) 
+         ENDIF
     
-         !                                                        ! module of wind stress and wind speed at T-point
-         zcoef = 1. / ( zrhoa * zcdrag )
-!CDIR NOVERRCHK
-         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
-               IF ( .NOT. (ln_shelf_flx .AND. ln_cpl)) THEN
-               wndm(ji,jj) = SQRT( zmod * zcoef )
-               ENDIF
-            END DO
-         END DO
          taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1)
-         CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
+         CALL lbc_lnk_multi( taum(:,:), 'T', 1., wndm(:,:), 'T', 1. )
 
          IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked)
