!!----------------------------------------------------------------------
   !!                      ***  ldfdyn_2d.h90  ***
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LODYC-IPSL  (2003)
   !!----------------------------------------------------------------------

   SUBROUTINE ldf_dyn_2d( ld_print )
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE ldf_dyn_2d  ***
      !!                  
      !! ** Purpose :   initializations of the horizontal ocean physics
      !!
      !! ** Method :
      !!             ORCA05 OCEAN VERSION with bilaplacian
      !!
      !!      2D eddy viscosity coefficients ( longitude, latitude )
      !!
      !!       biharmonic operator : ahm1 is defined at u-point
      !!                             ahm2 is defined at v-point
      !!                           : ahm3, ahm4 not used
      !!
      !!----------------------------------------------------------------------
      !! * Arguments
      LOGICAL, INTENT(in) :: ld_print   ! print coef. in numout

      !! * Local variables
      REAL(wp) ::   za00, zdx
      !!----------------------------------------------------------------------

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient for bilaplacian'
      IF(lwp) WRITE(numout,*) '~~~~~~  --'

      ! Read 2d integer array to specify western boundary increase in the
      ! ===================== equatorial strip (20N-20S) defined at t-points

      ! Set ahm1 and ahm2
      ! =================

      ! define ahm1 and ahm2 at the right grid point position
      ! biharmonic : ahm1 (ahm2) defined at u- (v-) point
      
      zdx = ra * rad * 0.5
      za00 = ahm0/(zdx*zdx*zdx)

      IF(lwp) WRITE(numout,*) '        orca05 ocean model eq resolution = ', zdx, ' eq ahm = ' ahm0
      IF(lwp) WRITE(numout,*) '        zdx=', zdx
      DO jj = 1, jpj
        DO ji = 1, jpi
          ahm1(ji,jj) = za00 * e1u(ji,jj) * e1u(ji,jj) * e1u(ji,jj)
          ahm2(ji,jj) = za00 * e1v(ji,jj) * e1v(ji,jj) * e1v(ji,jj)
        END DO
      END DO 

      ! Control print
      IF( lwp .AND. ld_print ) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'inildf: 2D ahm1 array'
         CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
         WRITE(numout,*)
         WRITE(numout,*) 'inildf: 2D ahm2 array'
         CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
      ENDIF

   END SUBROUTINE ldf_dyn_2d