 MODULE nudging
!!======================================================================
!!                   ***  MODULE  nudging  ***
!! Hydrodynamics: relaxation towards climatological, OGCM or zonal mean
!!======================================================================





























































! TO DO list
! - passer tous les logiques, timescale de relaxation etc dans paraspec.txt et
! cas_init_nudging pour cas test
! - ajouter l_nud_sal pour choisir quelle variable on nudge
!! * Modules used
   USE parameters, ONLY : rsh,rlg

   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC nud_init, &               ! called by main.F90
          nud_correct_fields, &     ! called by step.F90
          nud_dealloc, &
          nud_init_out, &           ! called by step.F90
          nud_out_cum               ! called by step.F90

!! * Shared or public module variables
   LOGICAL, PUBLIC       :: l_nud_read_init = .FALSE.
   LOGICAL, PUBLIC       :: l_nud_ssh = .TRUE.,      &
                            l_nud_temp = .TRUE.,      &
                            l_nud_sal = .FALSE.,      &
                            l_nud_boz = .FALSE.,      &
                            l_nud_u3d = .TRUE.,      &
                            l_nud_ubt = .TRUE.,      &
                            l_nud_zonalmean = .TRUE.
   LOGICAL, PUBLIC, ALLOCATABLE, DIMENSION(:) :: l_out_nud
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:,:)   :: ssh_clim, u_clim, v_clim
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:,:,:) :: temp_clim, sal_clim, boz_clim, &
                                                         uz_clim, vz_clim
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:)   :: ssh_clim_zonal, u_clim_zonal, v_clim_zonal, &
                                                       ssh_zonal, u_zonal, v_zonal
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:,:) :: temp_clim_zonal, sal_clim_zonal, boz_clim_zonal, &
                                                       uz_clim_zonal, vz_clim_zonal,uz_zonal, vz_zonal, &
                                                       temp_zonal, sal_zonal, boz_zonal
   REAL(KIND=rlg), PUBLIC   :: nud_dt_cum
   REAL(KIND=rsh)   :: nud_timescale

! Arrays to save averaged nudging field
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:,:)   :: ssh_clim_cum, u_clim_cum, v_clim_cum
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:,:,:) :: temp_clim_cum, sal_clim_cum, boz_clim_cum, &
                                                         uz_clim_cum, vz_clim_cum
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:)   :: ssh_zonal_cum, u_zonal_cum, v_zonal_cum
   REAL(KIND=rsh),PUBLIC,ALLOCATABLE,DIMENSION(:,:) :: uz_zonal_cum, vz_zonal_cum, &
                                                       temp_zonal_cum, sal_zonal_cum, boz_zonal_cum

   REAL(kind=rlg),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diagx_velnud,diagy_velnud
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:),PUBLIC   :: diag_sshnud

!! * Private variables
   INTEGER         :: comm_zonal

 CONTAINS


!!======================================================================

  SUBROUTINE nud_init(l_ssh,l_sal,l_temp,l_boz,l_u3d,l_ubt,l_zonalmean)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE nud_init  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2013-09  (V. Garnier) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : liminm1,liminm2,limaxp2,ljminm1,ljminm2,ljmaxp2,kmax, &
                          limin,limax
   USE comvarp2d, ONLY : sal,temp,ssh,u,v
   USE comvarp3d, ONLY : uz,vz,boz


!! * Arguments
   LOGICAL, INTENT( in ) :: l_ssh,l_sal,l_temp,l_boz,l_u3d,l_ubt,l_zonalmean

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter NUD_INIT routine'

   nud_timescale = 50.0_rsh*86400.0_rsh    ! TODO a lire

! Allocate climatological (or OGCM...) variables
   IF (l_ssh) ALLOCATE( ssh_clim(liminm2:limaxp2,ljminm2:ljmaxp2) )
   IF (l_ssh .AND. l_zonalmean) THEN
     ALLOCATE( ssh_clim_zonal(ljminm2:ljmaxp2),ssh_zonal(ljminm2:ljmaxp2) )
     ssh_clim_zonal(ljminm2:ljmaxp2)=ssh((limax+limin)/2,ljminm2:ljmaxp2)

     ssh_zonal(:)=0.0_rsh
   END IF
   IF (l_ubt) THEN
     ALLOCATE( u_clim(liminm2:limaxp2,ljminm2:ljmaxp2) )
     ALLOCATE( v_clim(liminm2:limaxp2,ljminm2:ljmaxp2) )
     IF (l_zonalmean) THEN
       ALLOCATE( u_clim_zonal(ljminm2:ljmaxp2),u_zonal(ljminm2:ljmaxp2) )
       ALLOCATE( v_clim_zonal(ljminm2:ljmaxp2),v_zonal(ljminm2:ljmaxp2) )
       u_clim_zonal(ljminm2:ljmaxp2)=u((limax+limin)/2,ljminm2:ljmaxp2)
       v_clim_zonal(ljminm2:ljmaxp2)=v((limax+limin)/2,ljminm2:ljmaxp2)

       u_zonal(:)=0.0_rsh
       v_zonal(:)=0.0_rsh
     END IF
   END IF
   IF (l_u3d) THEN
     ALLOCATE( uz_clim(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2) )
     ALLOCATE( vz_clim(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2) )
     IF (l_zonalmean) THEN
       ALLOCATE( uz_clim_zonal(1:kmax,ljminm2:ljmaxp2),uz_zonal(1:kmax,ljminm2:ljmaxp2) )
       ALLOCATE( vz_clim_zonal(1:kmax,ljminm2:ljmaxp2),vz_zonal(1:kmax,ljminm2:ljmaxp2) )
       uz_clim_zonal(1:kmax,ljminm2:ljmaxp2)=uz(1:kmax,(limax+limin)/2,ljminm2:ljmaxp2)
       vz_clim_zonal(1:kmax,ljminm2:ljmaxp2)=vz(1:kmax,(limax+limin)/2,ljminm2:ljmaxp2)

       uz_zonal(:,:)=0.0_rsh
       vz_zonal(:,:)=0.0_rsh
     END IF
   END IF
   IF (l_sal) ALLOCATE( sal_clim(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   IF (l_sal .AND. l_zonalmean) THEN
     ALLOCATE( sal_clim_zonal(1:kmax,ljminm1:ljmaxp2),sal_zonal(1:kmax,ljminm1:ljmaxp2) )
     sal_clim_zonal(1:kmax,ljminm1:ljmaxp2)=sal(1:kmax,(limax+limin)/2,ljminm1:ljmaxp2)
     sal_zonal(:,:)=0.0_rsh
   END IF
   IF (l_temp) ALLOCATE( temp_clim(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   IF (l_temp .AND. l_zonalmean) THEN
     ALLOCATE( temp_clim_zonal(1:kmax,ljminm1:ljmaxp2),temp_zonal(1:kmax,ljminm1:ljmaxp2) )
     temp_clim_zonal(1:kmax,ljminm1:ljmaxp2)=temp(1:kmax,(limax+limin)/2,ljminm1:ljmaxp2)

     temp_zonal(:,:)=0.0_rsh
   END IF
   IF (l_boz) ALLOCATE( boz_clim(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) )
   IF (l_boz .AND. l_zonalmean) THEN
     ALLOCATE( boz_clim_zonal(1:kmax,ljminm1:ljmaxp2),boz_zonal(1:kmax,ljminm1:ljmaxp2) )
     boz_clim_zonal(1:kmax,ljminm1:ljmaxp2)=boz(1:kmax,(limax+limin)/2,ljminm1:ljmaxp2)
     boz_zonal(:,:)=0.0_rsh
   END IF

   IF (l_nud_read_init) THEN
      CALL nud_init_read
   ELSE
      IF (l_ssh) ssh_clim(:,:)=ssh(:,:)
      IF (l_ubt) THEN
        u_clim(:,:)=u(:,:)
        v_clim(:,:)=v(:,:)
      END IF
      IF (l_sal) sal_clim(:,:,:)=sal(:,:,:)
      IF (l_temp) temp_clim(:,:,:)=temp(:,:,:)
      IF (l_boz) boz_clim(:,:,:)=boz(:,:,:)
      IF (l_u3d) THEN
        uz_clim(:,:,:)=uz(:,:,:)
        vz_clim(:,:,:)=vz(:,:,:)
      END IF
   END IF

   !PRINT*, 'exit NUD_INIT routine'

  END SUBROUTINE nud_init

!!======================================================================

  SUBROUTINE nud_dealloc

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE nud_dealloc  ***
!&E
!&E ** Purpose :  deallocate nudging variables
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2013-09  (V. Garnier) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Arguments

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter NUD_DEALLOC routine'

   IF (l_nud_ssh) DEALLOCATE( ssh_clim )
   IF (l_nud_ssh .AND. l_nud_zonalmean) DEALLOCATE( ssh_clim_zonal,ssh_zonal )
   IF (l_nud_ubt) THEN
     DEALLOCATE( u_clim )
     DEALLOCATE( v_clim )
     IF (l_nud_zonalmean) THEN
       DEALLOCATE( u_clim_zonal,u_zonal )
       DEALLOCATE( v_clim_zonal,v_zonal )
     END IF
   END IF
   IF (l_nud_u3d) THEN
     DEALLOCATE( uz_clim )
     DEALLOCATE( vz_clim )
     IF (l_nud_zonalmean) THEN
       DEALLOCATE( uz_clim_zonal,uz_zonal )
       DEALLOCATE( vz_clim_zonal,vz_zonal )
     END IF
   END IF
   IF (l_nud_sal) DEALLOCATE( sal_clim )
   IF (l_nud_sal .AND. l_nud_zonalmean) DEALLOCATE( sal_clim_zonal,sal_zonal )
   IF (l_nud_temp) DEALLOCATE( temp_clim )
   IF (l_nud_temp .AND. l_nud_zonalmean) DEALLOCATE( temp_clim_zonal,temp_zonal )
   IF (l_nud_boz) DEALLOCATE( boz_clim )
   IF (l_nud_boz .AND. l_nud_zonalmean) DEALLOCATE( boz_clim_zonal,boz_zonal )

   !PRINT*, 'exit NUD_DEALLOC routine'

  END SUBROUTINE nud_dealloc

!!======================================================================

  SUBROUTINE nud_init_read

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE nud_init_read  ***
!&E
!&E ** Purpose :   pour PIERRE lecture d un fichier
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2013-09  (V. Garnier) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Arguments

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter NUD_INIT_READ routine'

   !PRINT*, 'exit NUD_INIT_READ routine'

  END SUBROUTINE nud_init_read

!!======================================================================


!!======================================================================

  SUBROUTINE nud_correct_fields(l_ssh,l_sal,l_temp,l_boz,l_u3d,l_ubt,l_zonalmean)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE nud_correct_fields  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2013-09  (V. Garnier) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : liminm2,liminm1,limaxp2, &
                          ljminm2,ljminm1,ljmaxp2, &
                          liminp1,limaxm1, &
                          kmax,ljmin,ljmax
   USE comvars2d,  ONLY : halfdt,h0,h0fond,ig,id,nbouc
   USE comvarp2d,  ONLY : sal,temp,ssh,u,v
   USE comvarp3d,  ONLY : uz,vz,boz
   USE toolmpi,ONLY : CREATE_COMM_ZONAL,ADD_ZONAL_RSH,iproc,ex_i_rsh,ex_j_rsh
   USE parameters, ONLY : limin,limax

!! * Arguments
   LOGICAL, INTENT( in ) :: l_ssh,l_sal,l_temp,l_boz,l_u3d,l_ubt,l_zonalmean

!! * Local declarations
   INTEGER         :: i,j
   REAL(KIND=rsh)  :: nud_invtimescale

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter NUD_CORRECT_FIELDS routine'

   nud_invtimescale = 1.0_rsh/nud_timescale

!$OMP PARALLEL default(none) &
!$OMP& PRIVATE(i,j) &
!$OMP& SHARED(u,v,liminm2,liminm1,limaxp2,ljminm2,ljminm1,ljmaxp2, &
!$OMP&        liminp1,limaxm1,h0,ssh,h0fond,ig,id,nbouc, &
!$OMP&        l_ssh,l_sal,l_temp,l_boz,l_u3d,l_ubt,l_zonalmean,halfdt, &
!$OMP&        nud_invtimescale,ssh_clim,sal,sal_clim,temp,temp_clim, &
!$OMP&        boz,boz_clim,u_clim,v_clim,uz,uz_clim,vz,vz_clim, &
!$OMP&        ssh_zonal,ssh_clim_zonal,sal_zonal,sal_clim_zonal, &
!$OMP&        temp_zonal,temp_clim_zonal,boz_zonal,boz_clim_zonal, &
!$OMP&        u_zonal,u_clim_zonal,v_zonal,v_clim_zonal,uz_zonal, &


!$OMP&        comm_zonal, &


!$OMP&        uz_clim_zonal,vz_zonal,vz_clim_zonal,ljmin,ljmax)

   IF (l_zonalmean) THEN  ! relaxation towards a meridional section

! Define MPI communicator for zonal sum
     IF (nbouc==1) THEN
!$OMP master
       CALL CREATE_COMM_ZONAL(comm_zonal)
!$OMP end master
     END IF

     IF (l_ssh) THEN

!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         ssh_zonal(j)=0.0_rsh
         DO i=MAX0(liminp1,ig(j)),MIN0(limaxm1,id(j)-1)
           ssh_zonal(j)=ssh_zonal(j)+ssh(i,j)
         END DO
       END DO
!$OMP END DO

!$OMP barrier
!$OMP master
       CALL ADD_ZONAL_RSH(ssh_zonal,1,1,ljminm2,ljmaxp2,comm_zonal)

       ssh_zonal(:)=ssh_zonal(:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)

!$OMP end master
!$OMP barrier
!$OMP flush(ssh_zonal)

!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         DO i=liminm2,limaxp2
           IF (h0(i,j)+ssh(i,j) > h0fond) &
             ssh(i,j)= ssh(i,j) + halfdt*nud_invtimescale*(ssh_clim_zonal(j)-ssh_zonal(j))
         END DO
       END DO
!$OMP END DO


     END IF

     IF (l_sal) THEN

!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         sal_zonal(:,j)=0.0_rsh
         DO i=MAX0(liminp1,ig(j)),MIN0(limaxm1,id(j)-1)
           sal_zonal(:,j)=sal_zonal(:,j)+sal(:,i,j)
         END DO
       END DO
!$OMP END DO

!$OMP barrier
!$OMP master
       CALL ADD_ZONAL_RSH(sal_zonal,1,kmax,ljminm1,ljmaxp2,comm_zonal)

       sal_zonal(:,:)=sal_zonal(:,:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)

!$OMP end master
!$OMP barrier
!$OMP flush(sal_zonal)

!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         DO i=liminm1,limaxp2
           sal(:,i,j)= sal(:,i,j) + halfdt*nud_invtimescale*(sal_clim_zonal(:,j)-sal_zonal(:,j))
         END DO
       END DO
!$OMP END DO

     END IF

     IF (l_temp) THEN

!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         temp_zonal(:,j)=0.0_rsh
         DO i=MAX0(liminp1,ig(j)),MIN0(limaxm1,id(j)-1)
           temp_zonal(:,j)=temp_zonal(:,j)+temp(:,i,j)
         END DO
       END DO
!$OMP END DO

!$OMP barrier
!$OMP master
       CALL ADD_ZONAL_RSH(temp_zonal,1,kmax,ljminm1,ljmaxp2,comm_zonal)

       temp_zonal(:,:)=temp_zonal(:,:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)

!$OMP end master
!$OMP barrier
!$OMP flush(temp_zonal)

!$OMP DO SCHEDULE(static)
!print*,'j,temp(40,i,j),temp_clim_zonal(40,j),halfdt,nud_invtimescale'
       DO j=ljminm1,ljmaxp2
         DO i=liminm1,limaxp2
           temp(:,i,j)= temp(:,i,j) + halfdt*nud_invtimescale*(temp_clim_zonal(:,j)-temp_zonal(:,j))
!IF (i==12) print*,j,temp(40,i,j),temp_clim_zonal(40,j),temp_zonal(40,j),halfdt,nud_invtimescale
         END DO
       END DO
!$OMP END DO

     END IF

     IF (l_boz) THEN

!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         boz_zonal(:,j)=0.0_rsh
         DO i=MAX0(liminp1,ig(j)),MIN0(limaxm1,id(j)-1)
           boz_zonal(:,j)=boz_zonal(:,j)+boz(:,i,j)
         END DO
       END DO
!$OMP END DO

!$OMP barrier
!$OMP master
       CALL ADD_ZONAL_RSH(boz_zonal,1,kmax,ljminm1,ljmaxp2,comm_zonal)

       boz_zonal(:,:)=boz_zonal(:,:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)

!$OMP end master
!$OMP barrier
!$OMP flush(boz_zonal)

!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         DO i=liminm1,limaxp2
           boz(:,i,j)= boz(:,i,j) + halfdt*nud_invtimescale*(boz_clim_zonal(:,j)-boz_zonal(:,j))
         END DO
       END DO
!$OMP END DO

     END IF

     IF (l_ubt) THEN

!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         u_zonal(j)=0.0_rsh
         v_zonal(j)=0.0_rsh
         DO i=MAX0(liminp1,ig(j)),MIN0(limaxm1,id(j)-1)
           u_zonal(j)=u_zonal(j)+u(i,j)
           v_zonal(j)=v_zonal(j)+v(i,j)
         END DO
       END DO
!$OMP END DO

!$OMP barrier
!$OMP master
       CALL ADD_ZONAL_RSH(u_zonal,1,1,ljminm2,ljmaxp2,comm_zonal)
       CALL ADD_ZONAL_RSH(v_zonal,1,1,ljminm2,ljmaxp2,comm_zonal)

       u_zonal(:)=u_zonal(:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)
       v_zonal(:)=v_zonal(:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)

!$OMP end master
!$OMP barrier
!$OMP flush(u_zonal,v_zonal)

!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         DO i=liminm2,limaxp2
           u(i,j)= u(i,j) + halfdt*nud_invtimescale*(u_clim_zonal(j)-u_zonal(j))
           v(i,j)= v(i,j) + halfdt*nud_invtimescale*(v_clim_zonal(j)-v_zonal(j))
         END DO
       END DO
!$OMP END DO

     END IF

     IF (l_u3d) THEN

!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         uz_zonal(:,j)=0.0_rsh
         vz_zonal(:,j)=0.0_rsh
         DO i=MAX0(liminp1,ig(j)),MIN0(limaxm1,id(j)-1)
           uz_zonal(:,j)=uz_zonal(:,j)+uz(:,i,j)
           vz_zonal(:,j)=vz_zonal(:,j)+vz(:,i,j)
         END DO
       END DO
!$OMP END DO

!$OMP barrier
!$OMP master
       CALL ADD_ZONAL_RSH(uz_zonal,1,kmax,ljminm2,ljmaxp2,comm_zonal)
       CALL ADD_ZONAL_RSH(vz_zonal,1,kmax,ljminm2,ljmaxp2,comm_zonal)

       uz_zonal(:,:)=uz_zonal(:,:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)
       vz_zonal(:,:)=vz_zonal(:,:)/((id((ljmin+ljmax)/2)-1)-ig((ljmin+ljmax)/2)+1)

!$OMP end master
!$OMP barrier
!$OMP flush(uz_zonal,vz_zonal)

!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         DO i=liminm2,limaxp2
           uz(:,i,j)= uz(:,i,j) + halfdt*nud_invtimescale*(uz_clim_zonal(:,j)-uz_zonal(:,j))
           vz(:,i,j)= vz(:,i,j) + halfdt*nud_invtimescale*(vz_clim_zonal(:,j)-vz_zonal(:,j))
         END DO
       END DO
!$OMP END DO


     END IF

!!!!!$OMP master
!!!!CALL MPI_COMM_FREE(comm_zonal,ierr_mpi)
!!!!!$OMP end master

   ELSE  ! relaxation towards 3D fields

     IF (l_ssh) THEN
!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         DO i=liminm2,limaxp2
           IF (h0(i,j)+ssh(i,j) > h0fond) &
             ssh(i,j)= ssh(i,j) + halfdt*nud_invtimescale*(ssh_clim(i,j)-ssh(i,j))
         END DO
       END DO
!$OMP END DO
!$OMP barrier
!$OMP master
       CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
       CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(ssh)
     END IF

     IF (l_sal) THEN
!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         DO i=liminm1,limaxp2
           sal(:,i,j)= sal(:,i,j) + halfdt*nud_invtimescale*(sal_clim(:,i,j)-sal(:,i,j))
         END DO
       END DO
!$OMP END DO
!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(sal)
     END IF

     IF (l_temp) THEN
!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         DO i=liminm1,limaxp2
           temp(:,i,j)= temp(:,i,j) + halfdt*nud_invtimescale*(temp_clim(:,i,j)-temp(:,i,j))
         END DO
       END DO
!$OMP END DO
!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(temp)
     END IF

     IF (l_boz) THEN
!$OMP DO SCHEDULE(static)
       DO j=ljminm1,ljmaxp2
         DO i=liminm1,limaxp2
           boz(:,i,j)= boz(:,i,j) + halfdt*nud_invtimescale*(boz_clim(:,i,j)-boz(:,i,j))
         END DO
       END DO
!$OMP END DO
!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,boz(:,liminm1:limaxp2,ljminm1:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(boz)
     END IF

     IF (l_ubt) THEN
!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         DO i=liminm2,limaxp2
           u(i,j)= u(i,j) + halfdt*nud_invtimescale*(u_clim(i,j)-u(i,j))
           v(i,j)= v(i,j) + halfdt*nud_invtimescale*(v_clim(i,j)-v(i,j))
         END DO
       END DO
!$OMP END DO
!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(u,v)
     END IF

     IF (l_u3d) THEN
!$OMP DO SCHEDULE(static)
       DO j=ljminm2,ljmaxp2
         DO i=liminm2,limaxp2
           uz(:,i,j)= uz(:,i,j) + halfdt*nud_invtimescale*(uz_clim(:,i,j)-uz(:,i,j))
           vz(:,i,j)= vz(:,i,j) + halfdt*nud_invtimescale*(vz_clim(:,i,j)-vz(:,i,j))
         END DO
       END DO
!$OMP END DO
!$OMP barrier
!$OMP master
   CALL ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
!$OMP end master
!$OMP barrier
!$OMP flush(uz,vz)
     END IF

   END IF

!$OMP END PARALLEL

   !PRINT*, 'exit NUD_CORRECT_FIELDS routine'

  END SUBROUTINE nud_correct_fields


!!======================================================================

  SUBROUTINE nud_init_out(nb)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE nud_init_out  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2013-09  (V. Garnier) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters,ONLY : limin,limax,ljmin,ljmax,kmax

!! * Arguments
   INTEGER, INTENT(in) :: nb

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter NUD_INIT routine'

   nud_dt_cum=0.0_rlg

! Allocate climatological (or OGCM...) variables to cumulate data
   IF (l_nud_ssh .AND. l_nud_zonalmean) THEN
     ALLOCATE( ssh_zonal_cum(ljmin:ljmax) )
     ssh_zonal_cum(ljmin:ljmax)=0.0_rsh
   ELSE IF (l_nud_ssh) THEN
     ALLOCATE( ssh_clim_cum(limin:limax,ljmin:ljmax) )
     ssh_clim_cum(limin:limax,ljmin:ljmax)=0.0_rsh
   END IF
   IF (l_nud_ubt) THEN
     IF (l_nud_zonalmean) THEN
       ALLOCATE( u_zonal_cum(ljmin:ljmax) )
       ALLOCATE( v_zonal_cum(ljmin:ljmax) )
       u_zonal_cum(ljmin:ljmax)=0.0_rsh
       v_zonal_cum(ljmin:ljmax)=0.0_rsh
     ELSE
       ALLOCATE( u_clim_cum(limin:limax,ljmin:ljmax) )
       ALLOCATE( v_clim_cum(limin:limax,ljmin:ljmax) )
       u_clim_cum(limin:limax,ljmin:ljmax)=0.0_rsh
       v_clim_cum(limin:limax,ljmin:ljmax)=0.0_rsh
     END IF
   END IF
   IF (l_nud_u3d) THEN
     IF (l_nud_zonalmean) THEN
       ALLOCATE( uz_zonal_cum(1:kmax,ljmin:ljmax) )
       ALLOCATE( vz_zonal_cum(1:kmax,ljmin:ljmax) )
       uz_zonal_cum(1:kmax,ljmin:ljmax)=0.0_rsh
       vz_zonal_cum(1:kmax,ljmin:ljmax)=0.0_rsh
     ELSE
       ALLOCATE( uz_clim_cum(1:kmax,limin:limax,ljmin:ljmax) )
       ALLOCATE( vz_clim_cum(1:kmax,limin:limax,ljmin:ljmax) )
       uz_clim_cum(1:kmax,limin:limax,ljmin:ljmax)=0.0_rsh
       vz_clim_cum(1:kmax,limin:limax,ljmin:ljmax)=0.0_rsh
     END IF
   END IF
   IF (l_nud_sal .AND. l_nud_zonalmean) THEN
     ALLOCATE( sal_zonal_cum(1:kmax,ljmin:ljmax) )
     sal_zonal_cum(1:kmax,ljmin:ljmax)=0.0_rsh
   ELSE IF (l_nud_sal) THEN
     ALLOCATE( sal_clim_cum(1:kmax,limin:limax,ljmin:ljmax) )
     sal_clim_cum(1:kmax,limin:limax,ljmin:ljmax)=0.0_rsh
   END IF
   IF (l_nud_temp .AND. l_nud_zonalmean) THEN
     ALLOCATE( temp_zonal_cum(1:kmax,ljmin:ljmax) )
     temp_zonal_cum(1:kmax,ljmin:ljmax)=0.0_rsh
   ELSE IF (l_nud_temp) THEN
     ALLOCATE( temp_clim_cum(1:kmax,limin:limax,ljmin:ljmax) )
     temp_clim_cum(1:kmax,limin:limax,ljmin:ljmax)=0.0_rsh
   END IF
   IF (l_nud_boz .AND. l_nud_zonalmean) THEN
     ALLOCATE( boz_zonal_cum(1:kmax,ljmin:ljmax) )
     boz_zonal_cum(1:kmax,ljmin:ljmax)=0.0_rsh
   ELSE IF (l_nud_boz) THEN
     ALLOCATE( boz_clim_cum(1:kmax,limin:limax,ljmin:ljmax) )
     boz_clim_cum(1:kmax,limin:limax,ljmin:ljmax)=0.0_rsh
   END IF

! Create the output file saving the nudging terms
! -----------------------------------------------

   CALL out_nud_createfile(nb)

   !PRINT*, 'exit NUD_INIT_OUT routine'

  END SUBROUTINE nud_init_out
!!======================================================================

  SUBROUTINE nud_out_cum(nb,dtp,dtcum,dtnext)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE nud_out_cum  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2013-09  (V. Garnier) Original code
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : limin,limax,ljmin,ljmax,kmax
   USE comvars2d,  ONLY : nbouc


!! * Arguments
   INTEGER, INTENT(in) :: nb
   REAL(KIND=rlg), INTENT(IN) :: dtp
   REAL(KIND=rlg), INTENT(IN), OPTIONAL :: dtnext
   REAL(KIND=rlg), INTENT(INOUT), OPTIONAL :: dtcum

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter NUD_OUT_CUM routine'

! Accumulate the data in time
! ---------------------------

   IF (l_nud_ssh) THEN
     IF (l_nud_zonalmean) THEN
       ssh_zonal_cum(:) = ssh_zonal_cum(:) + dtp * ssh_zonal(ljmin:ljmax)
     ELSE 
       ssh_clim_cum(:,:)=ssh_clim_cum(:,:) + dtp * ssh_clim(limin:limax,ljmin:ljmax)
     END IF
   END IF

   IF (l_nud_sal) THEN
     IF (l_nud_zonalmean) THEN
       sal_zonal_cum(:,:) = sal_zonal_cum(:,:) + dtp * sal_zonal(1:kmax,ljmin:ljmax)
     ELSE 
       sal_clim_cum(:,:,:)= sal_clim_cum(:,:,:) + dtp * sal_clim(1:kmax,limin:limax,ljmin:ljmax)
     END IF
   END IF

   IF (l_nud_temp) THEN
     IF (l_nud_zonalmean) THEN
       temp_zonal_cum(:,:) = temp_zonal_cum(:,:) + dtp * temp_zonal(1:kmax,ljmin:ljmax)
     ELSE 
       temp_clim_cum(:,:,:)= temp_clim_cum(:,:,:) + dtp * temp_clim(1:kmax,limin:limax,ljmin:ljmax)
     END IF
   END IF

   IF (l_nud_boz) THEN
     IF (l_nud_zonalmean) THEN
       boz_zonal_cum(:,:) = boz_zonal_cum(:,:) + dtp * boz_zonal(1:kmax,ljmin:ljmax)
     ELSE 
       boz_clim_cum(:,:,:)= boz_clim_cum(:,:,:) + dtp * boz_clim(1:kmax,limin:limax,ljmin:ljmax)
     END IF
   END IF

   IF (l_nud_ubt) THEN
     IF (l_nud_zonalmean) THEN
       u_zonal_cum(:) = u_zonal_cum(:) + dtp * u_zonal(ljmin:ljmax)
       v_zonal_cum(:) = v_zonal_cum(:) + dtp * v_zonal(ljmin:ljmax)
     ELSE 
       u_clim_cum(:,:)= u_clim_cum(:,:) + dtp * u_clim(limin:limax,ljmin:ljmax)
       v_clim_cum(:,:)= v_clim_cum(:,:) + dtp * v_clim(limin:limax,ljmin:ljmax)
     END IF
   END IF

   IF (l_nud_u3d) THEN
     IF (l_nud_zonalmean) THEN
       uz_zonal_cum(:,:) = uz_zonal_cum(:,:) + dtp * uz_zonal(1:kmax,ljmin:ljmax)
       vz_zonal_cum(:,:) = vz_zonal_cum(:,:) + dtp * vz_zonal(1:kmax,ljmin:ljmax)
     ELSE 
       uz_clim_cum(:,:,:)= uz_clim_cum(:,:,:) + dtp * uz_clim(1:kmax,limin:limax,ljmin:ljmax)
       vz_clim_cum(:,:,:)= vz_clim_cum(:,:,:) + dtp * vz_clim(1:kmax,limin:limax,ljmin:ljmax)
     END IF
   END IF

   IF ( PRESENT(dtnext) .AND. nbouc>1 ) THEN

     IF (l_nud_ssh) THEN
       IF (l_nud_zonalmean) THEN
         ssh_zonal_cum(:) = ssh_zonal_cum(ljmin:ljmax) / dtcum
       ELSE 
         ssh_clim_cum(:,:)= ssh_clim_cum(limin:limax,ljmin:ljmax) / dtcum
       END IF
     END IF

     IF (l_nud_sal) THEN
       IF (l_nud_zonalmean) THEN
         sal_zonal_cum(:,:) = sal_zonal_cum(1:kmax,ljmin:ljmax) / dtcum
       ELSE 
         sal_clim_cum(:,:,:)= sal_clim_cum(1:kmax,limin:limax,ljmin:ljmax) / dtcum
       END IF
     END IF

     IF (l_nud_temp) THEN
       IF (l_nud_zonalmean) THEN
         temp_zonal_cum(:,:) = temp_zonal_cum(1:kmax,ljmin:ljmax) / dtcum
       ELSE 
         temp_clim_cum(:,:,:) = temp_clim_cum(1:kmax,limin:limax,ljmin:ljmax) / dtcum
       END IF
     END IF

     IF (l_nud_boz) THEN
       IF (l_nud_zonalmean) THEN
         boz_zonal_cum(:,:) = boz_zonal_cum(1:kmax,ljmin:ljmax) / dtcum
       ELSE 
         boz_clim_cum(:,:,:) = boz_clim_cum(1:kmax,limin:limax,ljmin:ljmax) / dtcum
       END IF
     END IF

     IF (l_nud_ubt) THEN
       IF (l_nud_zonalmean) THEN
         u_zonal_cum(:) = u_zonal_cum(ljmin:ljmax) / dtcum
         v_zonal_cum(:) = v_zonal_cum(ljmin:ljmax) / dtcum
       ELSE 
         u_clim_cum(:,:)= u_clim_cum(limin:limax,ljmin:ljmax) / dtcum
         v_clim_cum(:,:)= v_clim_cum(limin:limax,ljmin:ljmax) / dtcum
       END IF
     END IF

     IF (l_nud_u3d) THEN
       IF (l_nud_zonalmean) THEN
         uz_zonal_cum(:,:) = uz_zonal_cum(1:kmax,ljmin:ljmax) / dtcum
         vz_zonal_cum(:,:) = vz_zonal_cum(1:kmax,ljmin:ljmax) / dtcum
       ELSE 
         uz_clim_cum(:,:,:)= uz_clim_cum(1:kmax,limin:limax,ljmin:ljmax) / dtcum
         vz_clim_cum(:,:,:)= vz_clim_cum(1:kmax,limin:limax,ljmin:ljmax) / dtcum
       END IF
     END IF

! Save the fields
! ---------------

     CALL out_nud_save(nb)

! Reinitialize the arrays for the next period
! -------------------------------------------

     dtcum = 0.0_rlg

     IF (l_nud_ssh) THEN
       IF (l_nud_zonalmean) THEN
         ssh_zonal_cum(:) = dtnext * ssh_zonal(ljmin:ljmax)
       ELSE 
         ssh_clim_cum(:,:)= dtnext * ssh_clim(limin:limax,ljmin:ljmax)
       END IF
     END IF

     IF (l_nud_sal) THEN
       IF (l_nud_zonalmean) THEN
         sal_zonal_cum(:,:) = dtnext * sal_zonal(1:kmax,ljmin:ljmax)
       ELSE 
         sal_clim_cum(:,:,:)= dtnext * sal_clim(1:kmax,limin:limax,ljmin:ljmax)
       END IF
     END IF

     IF (l_nud_temp) THEN
       IF (l_nud_zonalmean) THEN
         temp_zonal_cum(:,:) =  dtnext * temp_zonal(1:kmax,ljmin:ljmax)
       ELSE 
         temp_clim_cum(:,:,:) = dtnext * temp_clim(1:kmax,limin:limax,ljmin:ljmax)
       END IF
     END IF

     IF (l_nud_boz) THEN
       IF (l_nud_zonalmean) THEN
         boz_zonal_cum(:,:) = dtnext * boz_zonal(1:kmax,ljmin:ljmax)
       ELSE 
         boz_clim_cum(:,:,:) = dtnext * boz_clim(1:kmax,limin:limax,ljmin:ljmax)
       END IF
     END IF

     IF (l_nud_ubt) THEN
       IF (l_nud_zonalmean) THEN
         u_zonal_cum(:) = dtnext * u_zonal(ljmin:ljmax)
         v_zonal_cum(:) = dtnext * v_zonal(ljmin:ljmax)
       ELSE 
         u_clim_cum(:,:)= dtnext * u_clim(limin:limax,ljmin:ljmax)
         v_clim_cum(:,:)= dtnext * v_clim(limin:limax,ljmin:ljmax)
       END IF
     END IF

     IF (l_nud_u3d) THEN
       IF (l_nud_zonalmean) THEN
         uz_zonal_cum(:,:) = dtnext * uz_zonal(1:kmax,ljmin:ljmax)
         vz_zonal_cum(:,:) = dtnext * vz_zonal(1:kmax,ljmin:ljmax)
       ELSE 
         uz_clim_cum(:,:,:)= dtnext * uz_clim(1:kmax,limin:limax,ljmin:ljmax)
         vz_clim_cum(:,:,:)= dtnext * vz_clim(1:kmax,limin:limax,ljmin:ljmax)
       END IF
     END IF

   END IF

   !PRINT*, 'exit NUD_OUT_CUM routine'

  END SUBROUTINE nud_out_cum

!!======================================================================

 
 END MODULE nudging
