Index: /NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaharm_fast.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaharm_fast.F90	(revision 15607)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaharm_fast.F90	(revision 15608)
@@ -83,4 +83,5 @@
    LOGICAL, PUBLIC :: ln_diaharm_fast
    LOGICAL, PUBLIC :: ln_diaharm_postproc_vel
+   LOGICAL, PUBLIC :: ln_diaharm_verbose
 
 
@@ -179,5 +180,5 @@
 
           
-          IF(lwp) WRITE(numout,*) 'diaharm_fast: sec2start = ',nint( (fjulday-fjulday_startharm)*86400._wp ),nsec_day - NINT(0.5_wp * rdt),adatrj * 86400._wp
+          IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) 'diaharm_fast: sec2start = ',nint( (fjulday-fjulday_startharm)*86400._wp ),nsec_day - NINT(0.5_wp * rdt),adatrj * 86400._wp
 
           IF( iom_use('tide_t') ) CALL iom_put( 'tide_t', sec2start )
@@ -366,6 +367,7 @@
       
 
-      NAMELIST/nam_diaharm_fast/ ln_diaharm_fast, ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, &
-               & tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily,ln_diaharm_postproc_vel
+      NAMELIST/nam_diaharm_fast/ ln_diaharm_fast, ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, &
+               & ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, &
+               & tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily,ln_diaharm_postproc_vel, ln_diaharm_verbose
       !!----------------------------------------------------------------------
       !JT
@@ -416,5 +418,6 @@
          WRITE(numout,*) '   Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily
          WRITE(numout,*) '   Number of Harmonics: nyear, nmonth = ', nyear, nmonth
-         WRITE(numout,*) '   Post-process velocity stats: ln_diaharm_postproc_vel = ', ln_diaharm_postproc_vel
+         WRITE(numout,*) '   Post-process velocity stats: ln_diaharm_postproc_vel = ',ln_diaharm_postproc_vel
+         WRITE(numout,*) '   Verbose: ln_diaharm_verbose = ', ln_diaharm_verbose
 
       ENDIF
@@ -447,7 +450,7 @@
       IF ( kt < 10 ) THEN
         ln_diaharm_read_restart = .FALSE.
-        IF(lwp) WRITE(numout,*) '   kt = ',kt
-        IF(lwp) WRITE(numout,*) '   kt < 10, so setting ln_diaharm_read_restart to .FALSE.'
-        IF(lwp) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart        
+        IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) '   kt = ',kt
+        IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) '   kt < 10, so setting ln_diaharm_read_restart to .FALSE.'
+        IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) '   Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart        
       ENDIF
 
@@ -777,8 +780,7 @@
 
       IF (ln_diaharm_postproc_vel)  THEN
-          IF (ln_ana_uvbar)  THEN
+          IF (ln_ana_uvbar) THEN
              ALLOCATE( amp_u2d(nb_ana,jpi,jpj), amp_v2d(nb_ana,jpi,jpj), phi_u2d(nb_ana,jpi,jpj), phi_v2d(nb_ana,jpi,jpj) )
-
-
+             
              ALLOCATE(tmp_u_amp_2d_mat(jpi,jpj),tmp_v_amp_2d_mat(jpi,jpj),tmp_u_phi_2d_mat(jpi,jpj),tmp_v_phi_2d_mat(jpi,jpj))
              ALLOCATE(a_u_2d_mat(jpi,jpj),b_u_2d_mat(jpi,jpj),a_v_2d_mat(jpi,jpj),b_v_2d_mat(jpi,jpj))
@@ -791,8 +793,7 @@
 
 
-          IF (ln_ana_uv3d)  THEN
+          IF (ln_ana_uv3d) THEN
              ALLOCATE( amp_u3d(nb_ana,jpi,jpj,jpk), amp_v3d(nb_ana,jpi,jpj,jpk), phi_u3d(nb_ana,jpi,jpj,jpk), phi_v3d(nb_ana,jpi,jpj,jpk) )
-
-
+             
              ALLOCATE(tmp_u_amp_3d_mat(jpi,jpj,jpk),tmp_v_amp_3d_mat(jpi,jpj,jpk),tmp_u_phi_3d_mat(jpi,jpj,jpk),tmp_v_phi_3d_mat(jpi,jpj,jpk))
              ALLOCATE(a_u_3d_mat(jpi,jpj,jpk),b_u_3d_mat(jpi,jpj,jpk),a_v_3d_mat(jpi,jpj,jpk),b_v_3d_mat(jpi,jpj,jpk))
@@ -842,17 +843,17 @@
              tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
              IF( iom_use(TRIM(tmp_name)) )  THEN
-                IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D)
-                IF(lwp) WRITE(numout,*) "diaharm_fast names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast names", tmp_name,tname(jh),' ',om_tide(jh), (2*rpi/3600.)/om_tide(jh),"hr"
                 CALL iom_put( TRIM(tmp_name), h_out2D(:,:) )
              ELSE
-                IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
              ENDIF
 
              tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
              IF( iom_use(TRIM(tmp_name)) )  THEN
-                IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)
                 CALL iom_put( TRIM(tmp_name), g_out2D(:,:) )
              ELSE
-                IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
              ENDIF
 
@@ -902,8 +903,8 @@
          tmp_name='TA_'//TRIM(suffix)//'_off'
          IF( iom_use(TRIM(tmp_name)) )  THEN
-            IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+            IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
             CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid))
          ELSE
-            IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
+            IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
          ENDIF
 
@@ -953,31 +954,17 @@
              tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
              IF( iom_use(TRIM(tmp_name)) )  THEN
-                IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)
                 CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) )
              ELSE
-                IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
              ENDIF
 
              tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
              IF( iom_use(TRIM(tmp_name)) )  THEN
-                IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)
                 CALL iom_put(tmp_name, g_out3D(:,:,:) )
              ELSE
-                IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
+                IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
              ENDIF
-
-          enddo             ! jh 
-
-         suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
-         tmp_name='TA_'//TRIM(suffix)//'_off'
-         IF( iom_use(TRIM(tmp_name)) )  THEN
-            IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
-            CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid))
-         ELSE
-            IF(lwp) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
-         ENDIF
-
-
-
 
 
@@ -1024,22 +1011,19 @@
              CALL FLUSH(numout)
 
+
+          enddo             ! jh 
+
+         suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
+         tmp_name='TA_'//TRIM(suffix)//'_off'
+         IF( iom_use(TRIM(tmp_name)) )  THEN
+            IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+            CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid))
+         ELSE
+            IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
+         ENDIF
+
       enddo                 ! jgrid
 
      CALL FLUSH(numout)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 
 
@@ -1078,25 +1062,19 @@
 
 
-                 DO jj = 2, nlcj !- 1
-                    DO ji = 2, nlci ! - 1
-
-    !             do jj=2,nlcj
-    !                do ji=2,nlci
-                        !IF ((ssumask(ji,jj) + ssumask(ji-1,jj)) == 0 ) CYCLE
-                        !IF ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) == 0 ) CYCLE
+                 DO jj = 1, nlcj !- 1
+                    DO ji = 1, nlci ! - 1
 
                         IF ( ((ssumask(ji,jj) + ssumask(ji-1,jj)) > 0 ) .AND. ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) > 0 ) ) THEN
-                            tmp_u_amp = ((amp_u2d(jh,ji,jj)*ssumask(ji,jj)) + (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj))
-                            tmp_v_amp = ((amp_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1))
-                            ! WORK ON THE WRAP AROUND
-                            !tmp_u_phi = ((phi_u2d(jh,ji,jj)*ssumask(ji,jj)) + (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj))
-                            !tmp_v_phi = ((phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1))
 
                             if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 1)) then
+
+                              tmp_u_amp = ((amp_u2d(jh,ji,jj)*ssumask(ji,jj)) + (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj))
                               !tmp_u_phi = ((phi_u2d(jh,ji,jj)*ssumask(ji,jj)) + (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj))
                               tmp_u_phi = atan2((sin(phi_u2d(jh,ji,jj)) + sin(phi_u2d(jh,ji-1,jj))),(cos(phi_u2d(jh,ji,jj)) + cos(phi_u2d(jh,ji-1,jj))))
                             else if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 0)) then
+                              tmp_u_amp = (amp_u2d(jh,ji,jj)*ssumask(ji,jj))
                               tmp_u_phi = (phi_u2d(jh,ji,jj)*ssumask(ji,jj))
                             else if ( (ssumask(ji,jj) == 0) .AND. (ssumask(ji-1,jj) == 1)) then
+                              tmp_u_amp = (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj))
                               tmp_u_phi = (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj))
                             else 
@@ -1106,95 +1084,18 @@
 
                             if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 1)) then
+                              tmp_v_amp = ((amp_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1))
                               !tmp_v_phi = ((phi_v2d(jh,ji,jj)*ssvmask(ji,jj)) + (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1))
                               tmp_v_phi = atan2((sin(phi_v2d(jh,ji,jj)) + sin(phi_v2d(jh,ji,jj-1))),(cos(phi_v2d(jh,ji,jj)) + cos(phi_v2d(jh,ji,jj-1))))
                             else if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 0)) then
+                              tmp_v_amp = (amp_v2d(jh,ji,jj)*ssvmask(ji,jj))
                               tmp_v_phi = (phi_v2d(jh,ji,jj)*ssvmask(ji,jj))
                             else if ( (ssvmask(ji,jj) == 0) .AND. (ssvmask(ji,jj-1) == 1)) then
-                              tmp_v_phi = (phi_v2d(ji,jj-1)*ssvmask(ji,jj-1))
+                              tmp_v_amp = (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1))
+                              !tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1))
+                              tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1))
                             else 
                               cycle
                             end if
 
-
-    !             do jj=1,nlcj
-    !                do ji=1,nlci
-
-    !                        tmp_u_amp = ((amp_u2d(jh,ji,jj)) + (amp_u2d(jh,ji-1,jj)))/(2.)
-    !                        tmp_v_amp = ((amp_v2d(jh,ji,jj)) + (amp_v2d(jh,ji,jj-1)))/(2.)
-    !                        ! WORK ON THE WRAP AROUND
-    !                        tmp_u_phi = ((phi_u2d(jh,ji,jj)) + (phi_u2d(jh,ji-1,jj)))/(2.)
-    !                        tmp_v_phi = ((phi_v2d(jh,ji,jj)) + (phi_v2d(jh,ji,jj-1)))/(2.)
-
-
-
-    !                        tmp_u_amp = (amp_u2d(jh,ji,jj)) 
-    !                        tmp_v_amp = (amp_v2d(jh,ji,jj)) 
-    !                        ! WORK ON THE WRAP AROUND
-    !                        tmp_u_phi = (phi_u2d(jh,ji,jj)) 
-    !                        tmp_v_phi = (phi_v2d(jh,ji,jj)) 
-
-
-
-!                            a_u = tmp_U_amp * cos(tmp_U_phi)
-!                            b_u = tmp_U_amp * sin(tmp_U_phi)
-!                            a_v = tmp_V_amp * cos(tmp_V_phi)
-!                            b_v = tmp_V_amp * sin(tmp_V_phi)
-
-!                            twodelta =  atan2( (tmp_V_amp**2  * sin( 2*(tmp_U_phi - tmp_V_phi)  ) ) , (   tmp_U_amp**2   +   tmp_V_amp**2  * cos( 2*(tmp_U_phi - tmp_V_phi)  )     ) )
-!                            delta = twodelta/2.
-
-!                            !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi))  )
-
-!                            tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 
-!                            if (tmpreal < 0) CYCLE
-!                            alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi))  )
-!                            if (alpha2 < 0) CYCLE
-!                            alpha= sqrt( alpha2 )
-
-
-!                            !major and minor axis of the ellipse
-!                            qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 )
-!                            !tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2
-!                            !qmin = 0
-!                            !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 )   ! but always positive.
-
-!                            tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2
-!                            if (tmpreal < 0) CYCLE
-!                            qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 )   ! but always positive.
-
-!                            !eccentricity of ellipse
-!                            tmpreal =  (qmax + qmin)
-!                            if (tmpreal < 0) CYCLE
-!                            ecc = (qmax - qmin)/(qmax + qmin)
-!                            ! Angle of major and minor ellipse
-!                            thetamax = atan2((  tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta)   ) , ( tmp_U_amp * cos( delta) )  )
-!                            thetamin = thetamax + rpi/2.
-
-
-
-!                            ! Rotary current components: Pugh A3.10
-!                            ! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors
-!                            ! so   Qc = clockwise     = anticyclonic = negative
-!                            ! and Qac = anticlockwise = cyclonic     = negative
-
-!                            tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
-!                            if (tmpreal < 0) CYCLE
-!                            Qc  = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))  )
-
-!                            tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 
-!                            if (tmpreal < 0) CYCLE
-!                            Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))  )
-
-
-!                            gc  = atan2(  (  (  tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  (  (tmp_U_amp*cos( tmp_U_phi ))  -  (tmp_V_amp*sin( tmp_V_phi ))  )  )
-!                            gac = atan2(  (  ( -tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  (  (tmp_U_amp*cos( tmp_U_phi ))  +  (tmp_V_amp*sin( tmp_V_phi ))  )  )
-
-!                            !Pugh A3.2
-!                            Phi_Ua = -0.5*(gac - gc)
-!                            dir_Ua = 0.5*(gac + gc)  ! positive from x axis
-
-!                            tmpreal = qmax
-!                            if (tmpreal < 0) CYCLE
-!                            polarity = (Qac - Qc)/qmax
 
                             a_u = tmp_U_amp * cos(tmp_U_phi)
@@ -1207,7 +1108,7 @@
 
                             !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi))  )
-
                             tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 
                             if (tmpreal < 0) tmpreal = 0 !CYCLE
+
                             !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi))  )
                             alpha2 = sqrt( tmpreal )
@@ -1218,7 +1119,4 @@
                             !major and minor axis of the ellipse
                             qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 )
-                            !tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2
-                            !qmin = 0
-                            !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 )   ! but always positive.
 
                             tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2
@@ -1228,9 +1126,9 @@
 
                             !eccentricity of ellipse
-
                             tmpreal =  (qmax + qmin)
-                            if (tmpreal == 0) tmpreal = tmpreal + 0.0000001! CYCLE
+                            if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE
                             !ecc = (qmax - qmin)/(qmax + qmin)
                             ecc = (qmax - qmin)/(tmpreal)
+
                             ! Angle of major and minor ellipse
                             thetamax = atan2((  tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta)   ) , ( tmp_U_amp * cos( delta) )  )
@@ -1245,16 +1143,18 @@
 
                             tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
-                            if (tmpreal < 0) tmpreal = 0! CYCLE
+                            if (tmpreal < 0) tmpreal = 0 !CYCLE
                             !Qc  = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))  )
                             Qc  = 0.5*sqrt( tmpreal )
 
                             tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 
-                            if (tmpreal < 0) tmpreal = 0! CYCLE
+                            if (tmpreal < 0) tmpreal = 0 !CYCLE
                             !Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))  )
                             Qac = 0.5*sqrt( tmpreal )
 
 
-                            gc  = atan2(  (  (  tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  (  (tmp_U_amp*cos( tmp_U_phi ))  -  (tmp_V_amp*sin( tmp_V_phi ))  )  )
-                            gac = atan2(  (  ( -tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  (  (tmp_U_amp*cos( tmp_U_phi ))  +  (tmp_V_amp*sin( tmp_V_phi ))  )  )
+                            gc  = atan2(  (  (  tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  &
+                                          & (  (tmp_U_amp*cos( tmp_U_phi ))  -  (tmp_V_amp*sin( tmp_V_phi ))  )  )
+                            gac = atan2(  (  ( -tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  &
+                                          & (  (tmp_U_amp*cos( tmp_U_phi ))  +  (tmp_V_amp*sin( tmp_V_phi ))  )  )
 
                             !Pugh A3.2
@@ -1270,5 +1170,4 @@
                             endif
 
-
                             tmp_u_amp_2d_mat(ji,jj) = tmp_u_amp
                             tmp_v_amp_2d_mat(ji,jj) = tmp_v_amp
@@ -1305,20 +1204,20 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                   IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                   IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                    CALL iom_put( TRIM(tmp_name), tmp_u_amp_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), tmp_v_amp_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), tmp_u_phi_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), tmp_v_phi_2d_mat(:,:))
                 ENDIF
@@ -1328,20 +1227,20 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                   IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                   IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                    CALL iom_put( TRIM(tmp_name), a_u_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), a_v_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), b_u_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), b_v_2d_mat(:,:))
                 ENDIF
@@ -1349,10 +1248,10 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), qmax_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), qmin_2d_mat(:,:))
                 ENDIF
@@ -1360,15 +1259,15 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), ecc_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), thetamax_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), thetamin_2d_mat(:,:))
                 ENDIF
@@ -1376,20 +1275,20 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), Qc_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), Qac_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), gc_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), gac_2d_mat(:,:))
                 ENDIF
@@ -1398,15 +1297,15 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), Phi_Ua_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), dir_Ua_2d_mat(:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), polarity_2d_mat(:,:))
                 ENDIF
@@ -1478,10 +1377,14 @@
 
 
-                  DO jk=1,jpkm1
+                DO jk=1,jpkm1
                      !DO jj = 2, nlcj ! - 1
                      !   DO ji = 2, nlci ! - 1
 
-                     DO jj = 2, jpjm1
-                        DO ji = 2, jpim1
+                 !    DO jj = 1, jpjm1    #works
+                 !       DO ji = 1, jpim1
+
+                     DO jj = 1, nlcj !- 1
+                        DO ji = 1, nlci ! - 1
+
 
         !             do jj=2,nlcj
@@ -1491,6 +1394,6 @@
 
                             IF ( ((umask(ji,jj,jk) + umask(ji-1,jj,jk)) > 0 ) .AND. ((vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) > 0 ) ) THEN
-                                tmp_u_amp = ((amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
-                                tmp_v_amp = ((amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
+                                !tmp_u_amp = ((amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
+                                !tmp_v_amp = ((amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
                                 ! WORK ON THE WRAP AROUND
                                 !tmp_u_phi = ((phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
@@ -1498,43 +1401,31 @@
 
                                 if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 1)) then
+                                  tmp_u_amp = ((amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
                                   !tmp_u_phi = ((phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) + (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
                                   tmp_u_phi = atan2((sin(phi_u3d(jh,ji,jj,jk)) + sin(phi_u3d(jh,ji-1,jj,jk))),(cos(phi_u3d(jh,ji,jj,jk)) + cos(phi_u3d(jh,ji-1,jj,jk))))
                                 else if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 0)) then
+                                  tmp_u_amp = (amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk)) 
                                   tmp_u_phi = (phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk))
                                 else if ( (umask(ji,jj,jk) == 0) .AND. (umask(ji-1,jj,jk) == 1)) then
+                                  tmp_u_amp = (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk))
                                   tmp_u_phi = (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk))
                                 else 
                                   cycle
-                                end if
+                                end if  
 
 
                                 if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 1)) then
+                                  tmp_v_amp = ((amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
                                   !tmp_v_phi = ((phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk)) + (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
                                   tmp_v_phi = atan2((sin(phi_v3d(jh,ji,jj,jk)) + sin(phi_v3d(jh,ji,jj-1,jk))),(cos(phi_v3d(jh,ji,jj,jk)) + cos(phi_v3d(jh,ji,jj-1,jk))))
                                 else if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 0)) then
+                                  tmp_v_amp = (amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk))
                                   tmp_v_phi = (phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk))
                                 else if ( (vmask(ji,jj,jk) == 0) .AND. (vmask(ji,jj-1,jk) == 1)) then
-                                  tmp_v_phi = (phi_v3d(ji,jj-1,jk)*vmask(ji,jj-1,jk))
+                                  tmp_v_amp = (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk))
+                                  tmp_v_phi = (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk))
                                 else 
                                   cycle
                                 end if
-
-
-        !             do jj=1,nlcj
-        !                do ji=1,nlci
-
-        !                        tmp_u_amp = ((amp_u2d(jh,ji,jj)) + (amp_u2d(jh,ji-1,jj)))/(2.)
-        !                        tmp_v_amp = ((amp_v2d(jh,ji,jj)) + (amp_v2d(jh,ji,jj-1)))/(2.)
-        !                        ! WORK ON THE WRAP AROUND
-        !                        tmp_u_phi = ((phi_u2d(jh,ji,jj)) + (phi_u2d(jh,ji-1,jj)))/(2.)
-        !                        tmp_v_phi = ((phi_v2d(jh,ji,jj)) + (phi_v2d(jh,ji,jj-1)))/(2.)
-
-
-
-        !                        tmp_u_amp = (amp_u2d(jh,ji,jj)) 
-        !                        tmp_v_amp = (amp_v2d(jh,ji,jj)) 
-        !                        ! WORK ON THE WRAP AROUND
-        !                        tmp_u_phi = (phi_u2d(jh,ji,jj)) 
-        !                        tmp_v_phi = (phi_v2d(jh,ji,jj)) 
 
 
@@ -1549,7 +1440,7 @@
 
                                 !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi))  )
-
                                 tmpreal = tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi)) 
                                 if (tmpreal < 0) tmpreal = 0 !CYCLE
+
                                 !alpha2 = sqrt( tmp_U_amp**4 + tmp_V_amp**4 + 2*tmp_U_amp**2*tmp_V_amp**2*cos(2*(tmp_U_phi - tmp_V_phi))  )
                                 alpha2 = sqrt( tmpreal )
@@ -1560,7 +1451,4 @@
                                 !major and minor axis of the ellipse
                                 qmax = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 + alpha**2)/2 )
-                                !tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2
-                                !qmin = 0
-                                !if (tmpreal > 0) qmin = sqrt( (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2 )   ! but always positive.
 
                                 tmpreal =  (tmp_U_amp**2 + tmp_V_amp**2 - alpha**2)/2
@@ -1570,9 +1458,9 @@
 
                                 !eccentricity of ellipse
-
                                 tmpreal =  (qmax + qmin)
-                                if (tmpreal == 0) tmpreal = tmpreal + 0.0000001! CYCLE
+                                if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE
                                 !ecc = (qmax - qmin)/(qmax + qmin)
                                 ecc = (qmax - qmin)/(tmpreal)
+
                                 ! Angle of major and minor ellipse
                                 thetamax = atan2((  tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta)   ) , ( tmp_U_amp * cos( delta) )  )
@@ -1587,16 +1475,18 @@
 
                                 tmpreal = tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
-                                if (tmpreal < 0) tmpreal = 0! CYCLE
+                                if (tmpreal < 0) tmpreal = 0 !CYCLE
                                 !Qc  = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 - (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))  )
                                 Qc  = 0.5*sqrt( tmpreal )
 
                                 tmpreal = tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi)) 
-                                if (tmpreal < 0) tmpreal = 0! CYCLE
+                                if (tmpreal < 0) tmpreal = 0 !CYCLE
                                 !Qac = 0.5*sqrt( tmp_U_amp**2 + tmp_V_amp**2 + (2*tmp_U_amp*tmp_V_amp*sin( tmp_V_phi - tmp_U_phi))  )
                                 Qac = 0.5*sqrt( tmpreal )
 
 
-                                gc  = atan2(  (  (  tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  (  (tmp_U_amp*cos( tmp_U_phi ))  -  (tmp_V_amp*sin( tmp_V_phi ))  )  )
-                                gac = atan2(  (  ( -tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  (  (tmp_U_amp*cos( tmp_U_phi ))  +  (tmp_V_amp*sin( tmp_V_phi ))  )  )
+                                gc  = atan2(  (  (  tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  &
+                                              & (  (tmp_U_amp*cos( tmp_U_phi ))  -  (tmp_V_amp*sin( tmp_V_phi ))  )  )
+                                gac = atan2(  (  ( -tmp_U_amp*sin( tmp_U_phi ) ) +  (tmp_V_amp*cos( tmp_V_phi)  ) )  ,  &
+                                              & (  (tmp_U_amp*cos( tmp_U_phi ))  +  (tmp_V_amp*sin( tmp_V_phi ))  )  )
 
                                 !Pugh A3.2
@@ -1613,5 +1503,4 @@
 
 
-
                                 tmp_u_amp_3d_mat(ji,jj,jk) = tmp_u_amp
                                 tmp_v_amp_3d_mat(ji,jj,jk) = tmp_v_amp
@@ -1642,27 +1531,27 @@
 
                             ENDIF
-                        END DO
-                     END DO
-                 END DO
+                        END DO !ji
+                     END DO    !jj
+                 END DO        !jk
 
 
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                   IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                   IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                    CALL iom_put( TRIM(tmp_name), tmp_u_amp_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), tmp_v_amp_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), tmp_u_phi_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), tmp_v_phi_3d_mat(:,:,:))
                 ENDIF
@@ -1672,20 +1561,20 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                   IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                   IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                    CALL iom_put( TRIM(tmp_name), a_u_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), a_v_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), b_u_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), b_v_3d_mat(:,:,:))
                 ENDIF
@@ -1693,10 +1582,10 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), qmax_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), qmin_3d_mat(:,:,:))
                 ENDIF
@@ -1704,15 +1593,15 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), ecc_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), thetamax_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), thetamin_3d_mat(:,:,:))
                 ENDIF
@@ -1720,20 +1609,20 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), Qc_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), Qac_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), gc_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), gac_3d_mat(:,:,:))
                 ENDIF
@@ -1742,15 +1631,15 @@
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), Phi_Ua_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), dir_Ua_3d_mat(:,:,:))
                 ENDIF
                 tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uv3d'
                 IF( iom_use(TRIM(tmp_name)) ) THEN
-                  IF(lwp) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
+                  IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
                   CALL iom_put( TRIM(tmp_name), polarity_3d_mat(:,:,:))
                 ENDIF
@@ -1770,5 +1659,5 @@
 
                 ecc_3d_mat(:,:,:) = 0.
-                thetamax_3d_mat(:,:,:) =0.
+                thetamax_3d_mat(:,:,:) = 0.
                 thetamin_3d_mat(:,:,:) = 0.
 
@@ -1783,5 +1672,5 @@
 
 
-             END DO
+             END DO    !jh
 
 
