Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 4682)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 4683)
@@ -18,4 +18,5 @@
    USE daymod
    USE tide_mod
+   !
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
@@ -34,9 +35,9 @@
    INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24
 
-   !                            !!!namelist variables
+   !                         !!** namelist variables **
    INTEGER ::   nit000_han    ! First time step used for harmonic analysis
    INTEGER ::   nitend_han    ! Last time step used for harmonic analysis
    INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis
-   INTEGER ::   nb_ana           ! Number of harmonics to analyse
+   INTEGER ::   nb_ana        ! Number of harmonics to analyse
 
    INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
@@ -119,5 +120,5 @@
             ENDIF
          END DO
-      ENDDO
+      END DO
       !
       IF(lwp) THEN
@@ -158,5 +159,5 @@
       ! ----------------------------
       ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
-      ana_temp(:,:,:,:) = 0.e0
+      ana_temp(:,:,:,:) = 0._wp
 
    END SUBROUTINE dia_harm_init
@@ -179,32 +180,31 @@
       IF( nn_timing == 1 )   CALL timing_start('dia_harm')
 
-      IF ( kt == nit000 ) CALL dia_harm_init
-
-      IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. &
-           (MOD(kt,nstep_han).EQ.0) ) THEN
-
-        ztime = (kt-nit000+1)*rdt 
+      IF( kt == nit000 ) CALL dia_harm_init
+
+      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
+
+         ztime = (kt-nit000+1) * rdt 
        
-        nhc = 0
-        DO jh = 1,nb_ana
-          DO jc = 1,2
-            nhc = nhc+1
-            ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
-                    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
-
-            DO jj = 1,jpj
-              DO ji = 1,jpi
-                ! Elevation
-                ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)        
+         nhc = 0
+         DO jh = 1, nb_ana
+            DO jc = 1, 2
+               nhc = nhc+1
+               ztemp =(     MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh))  &
+                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))
+
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ! Elevation
+                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask(ji,jj,1)
 #if defined key_dynspg_ts
-                ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
-                ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
-#endif
-              END DO
-            END DO
-
-          END DO
-        END DO
-       
+                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)
+                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)
+#endif
+                  END DO
+               END DO
+               !
+            END DO
+         END DO
+         !       
       END IF
 
@@ -249,6 +249,6 @@
          keq = keq + 1
          kun = 0
-         DO jh = 1,nb_ana
-            DO jc = 1,2
+         DO jh = 1, nb_ana
+            DO jc = 1, 2
                kun = kun + 1
                ksp = ksp + 1
@@ -296,7 +296,7 @@
                out_eta(ji,jj,jh       ) = X1 * tmask(ji,jj,1)
                out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1)
-            ENDDO
-         ENDDO
-      ENDDO
+            END DO
+         END DO
+      END DO
 
       ! ubar:
@@ -309,6 +309,6 @@
                   kun = kun + 1
                   ztmp4(kun)=ana_temp(ji,jj,kun,2)
-               ENDDO
-            ENDDO
+               END DO
+            END DO
 
             CALL SUR_DETERMINE(jj+1)
@@ -316,6 +316,6 @@
             ! Fill output array
             DO jh = 1, nb_ana
-               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
-               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
+               ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
+               ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
             END DO
 
@@ -326,11 +326,11 @@
          DO ji = 1, jpi
             DO jh = 1, nb_ana 
-               X1=ana_amp(ji,jj,jh,1)
-               X2=-ana_amp(ji,jj,jh,2)
-               out_u(ji,jj,jh) = X1 * umask(ji,jj,1)
-               out_u (ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
-            ENDDO
-         ENDDO
-      ENDDO
+               X1 = ana_amp(ji,jj,jh,1)
+               X2 =-ana_amp(ji,jj,jh,2)
+               out_u(ji,jj,jh       ) = X1 * umask(ji,jj,1)
+               out_u(ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1)
+            END DO
+         END DO
+      END DO
 
       ! vbar:
@@ -343,6 +343,6 @@
                   kun = kun + 1
                   ztmp4(kun)=ana_temp(ji,jj,kun,3)
-               ENDDO
-            ENDDO
+               END DO
+            END DO
 
             CALL SUR_DETERMINE(jj+1)
@@ -364,7 +364,7 @@
                out_v(ji,jj,jh)=X1 * vmask(ji,jj,1)
                out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1)
-            ENDDO
-         ENDDO
-      ENDDO
+            END DO
+         END DO
+      END DO
 
       CALL dia_wri_harm ! Write results in files
@@ -437,118 +437,118 @@
 #else
       DO jh = 1, nb_ana
-         CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh       ) )
-         CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,jh+nb_ana) )
-      END DO
-#endif
-
+         CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) )
+         CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
+      END DO
+#endif
+      !
    END SUBROUTINE dia_wri_harm
 
 
    SUBROUTINE SUR_DETERMINE(init)
-   !!---------------------------------------------------------------------------------
-   !!                      *** ROUTINE SUR_DETERMINE ***
-   !!    
-   !!    
-   !!       
-   !!---------------------------------------------------------------------------------
-   INTEGER, INTENT(in) ::   init 
-   !
-   INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
-   REAL(wp)                        :: zval1, zval2, zx1
-   REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
-   INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
-   !---------------------------------------------------------------------------------
-   CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
-   CALL wrk_alloc( jpincomax , ipos2 , ipivot        )
+      !!---------------------------------------------------------------------------------
+      !!                      *** ROUTINE SUR_DETERMINE ***
+      !!    
+      !!    
+      !!       
+      !!---------------------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   init 
+      !
+      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
+      REAL(wp)                        :: zval1, zval2, zx1
+      REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
+      INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
+      !---------------------------------------------------------------------------------
+      CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
+      CALL wrk_alloc( jpincomax , ipos2 , ipivot        )
             
-   IF( init == 1 ) THEN
-      IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
-      IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
-      !
-      ztmp3(:,:) = 0._wp
-      !
-      DO jk1_sd = 1, nsparse
-         DO jk2_sd = 1, nsparse
-            nisparse(jk2_sd) = nisparse(jk2_sd)
-            njsparse(jk2_sd) = njsparse(jk2_sd)
-            IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
-               ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
-                                                        + valuesparse(jk1_sd)*valuesparse(jk2_sd)
-            ENDIF
-         END DO
-      END DO
-
-      DO jj_sd = 1 ,ninco
-          ipos1(jj_sd) = jj_sd
-          ipos2(jj_sd) = jj_sd
-      ENDDO
-
-      DO ji_sd = 1 , ninco
-
-         !find greatest non-zero pivot:
-         zval1 = ABS(ztmp3(ji_sd,ji_sd))
-
-         ipivot(ji_sd) = ji_sd
-         DO jj_sd = ji_sd, ninco
-            zval2 = ABS(ztmp3(ji_sd,jj_sd))
-            IF( zval2.GE.zval1 )THEN
-               ipivot(ji_sd) = jj_sd
-               zval1         = zval2
-            ENDIF
-         ENDDO
-
-         DO ji1_sd = 1, ninco
-            zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
-            zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
-            ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
-            ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
-         ENDDO
-
-         ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
-         ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
-         ipos1(ji_sd)         = ipos2(ji_sd)
-         ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
-         zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
-         DO jj_sd = 1, ninco
-            ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
-         ENDDO
-
+      IF( init == 1 ) THEN
+         IF( nsparse > jpdimsparse )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
+         IF( ninco   > jpincomax   )   CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
+         !
+         ztmp3(:,:) = 0._wp
+         !
+         DO jk1_sd = 1, nsparse
+            DO jk2_sd = 1, nsparse
+               nisparse(jk2_sd) = nisparse(jk2_sd)
+               njsparse(jk2_sd) = njsparse(jk2_sd)
+               IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
+                  ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
+                     &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd)
+               ENDIF
+            END DO
+         END DO
+         !
+         DO jj_sd = 1 ,ninco
+            ipos1(jj_sd) = jj_sd
+            ipos2(jj_sd) = jj_sd
+         END DO
+         !
+         DO ji_sd = 1 , ninco
+            !
+            !find greatest non-zero pivot:
+            zval1 = ABS(ztmp3(ji_sd,ji_sd))
+            !
+            ipivot(ji_sd) = ji_sd
+            DO jj_sd = ji_sd, ninco
+               zval2 = ABS(ztmp3(ji_sd,jj_sd))
+               IF( zval2.GE.zval1 )THEN
+                  ipivot(ji_sd) = jj_sd
+                  zval1         = zval2
+               ENDIF
+            END DO
+            !
+            DO ji1_sd = 1, ninco
+               zcol1(ji1_sd)               = ztmp3(ji1_sd,ji_sd)
+               zcol2(ji1_sd)               = ztmp3(ji1_sd,ipivot(ji_sd))
+               ztmp3(ji1_sd,ji_sd)         = zcol2(ji1_sd)
+               ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
+            END DO
+            !
+            ipos2(ji_sd)         = ipos1(ipivot(ji_sd))
+            ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
+            ipos1(ji_sd)         = ipos2(ji_sd)
+            ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
+            zpivot(ji_sd)        = ztmp3(ji_sd,ji_sd)
+            DO jj_sd = 1, ninco
+               ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
+            END DO
+            !
+            DO ji2_sd = ji_sd+1, ninco
+               zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
+               DO jj_sd=1,ninco
+                  ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
+               END DO
+            END DO
+            !
+         END DO
+         !
+      ENDIF ! End init==1
+
+      DO ji_sd = 1, ninco
+         ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
          DO ji2_sd = ji_sd+1, ninco
-            zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
-            DO jj_sd=1,ninco
-               ztmp3(ji2_sd,jj_sd)=  ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
-            ENDDO
-         ENDDO
-
-      ENDDO
-
-   ENDIF ! End init==1
-
-   DO ji_sd = 1, ninco
-      ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
-      DO ji2_sd = ji_sd+1, ninco
-         ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
-      ENDDO
-   ENDDO
-
-   !system solving: 
-   ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
-   ji_sd = ninco
-   DO ji_sd = ninco-1, 1, -1
-      zx1=0.
-      DO jj_sd = ji_sd+1, ninco
-         zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
-      ENDDO
-      ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
-   ENDDO
-
-   DO jj_sd =1, ninco
-      ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
-   ENDDO
-
-   CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
-   CALL wrk_dealloc( jpincomax , ipos2 , ipivot        )
-
-  END SUBROUTINE SUR_DETERMINE
+            ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
+         END DO
+      END DO
+
+      !system solving: 
+      ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
+      ji_sd = ninco
+      DO ji_sd = ninco-1, 1, -1
+         zx1 = 0._wp
+         DO jj_sd = ji_sd+1, ninco
+            zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
+         END DO
+         ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
+      END DO
+
+      DO jj_sd =1, ninco
+         ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
+      END DO
+
+      CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
+      CALL wrk_dealloc( jpincomax , ipos2 , ipivot        )
+      !
+   END SUBROUTINE SUR_DETERMINE
 
 #else
