Index: /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90	(revision 13180)
+++ /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90	(revision 13181)
@@ -442,5 +442,5 @@
                
                zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0)
-                  &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts
+                  &            - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp  * ztmelts
 
                zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0)
Index: /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icevar.F90	(revision 13180)
+++ /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icevar.F90	(revision 13181)
@@ -635,14 +635,15 @@
       !!-------------------------------------------------------------------
       !
-      WHERE( pa_i (1:npti,:)   < 0._wp .AND. pa_i (1:npti,:)   > -epsi10 )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
-      WHERE( pv_i (1:npti,:)   < 0._wp .AND. pv_i (1:npti,:)   > -epsi10 )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
-      WHERE( pv_s (1:npti,:)   < 0._wp .AND. pv_s (1:npti,:)   > -epsi10 )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
-      WHERE( psv_i(1:npti,:)   < 0._wp .AND. psv_i(1:npti,:)   > -epsi10 )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
-      WHERE( poa_i(1:npti,:)   < 0._wp .AND. poa_i(1:npti,:)   > -epsi10 )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
-      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
-      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
+
+      WHERE( pa_i (1:npti,:)   < 0._wp )   pa_i (1:npti,:)   = 0._wp   !  a_i must be >= 0
+      WHERE( pv_i (1:npti,:)   < 0._wp )   pv_i (1:npti,:)   = 0._wp   !  v_i must be >= 0
+      WHERE( pv_s (1:npti,:)   < 0._wp )   pv_s (1:npti,:)   = 0._wp   !  v_s must be >= 0
+      WHERE( psv_i(1:npti,:)   < 0._wp )   psv_i(1:npti,:)   = 0._wp   ! sv_i must be >= 0
+      WHERE( poa_i(1:npti,:)   < 0._wp )   poa_i(1:npti,:)   = 0._wp   ! oa_i must be >= 0
+      WHERE( pe_i (1:npti,:,:) < 0._wp )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0
+      WHERE( pe_s (1:npti,:,:) < 0._wp )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0
       IF( ln_pnd_H12 ) THEN
-         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
-         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+         WHERE( pa_ip(1:npti,:) < 0._wp )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
+         WHERE( pv_ip(1:npti,:) < 0._wp )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
       ENDIF
       !
Index: /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/ICB/icblbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/ICB/icblbc.F90	(revision 13180)
+++ /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/ICB/icblbc.F90	(revision 13181)
@@ -81,5 +81,4 @@
       TYPE(iceberg), POINTER ::   this
       TYPE(point)  , POINTER ::   pt
-      INTEGER                ::   iine
       !!----------------------------------------------------------------------
 
@@ -92,8 +91,7 @@
          DO WHILE( ASSOCIATED(this) )
             pt => this%current_point
-            iine = INT( pt%xi + 0.5 )
-            IF( iine > mig(nicbei) ) THEN
+            IF( pt%xi > REAL(mig(nicbei),wp) + 0.5_wp ) THEN
                pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp
-            ELSE IF( iine < mig(nicbdi) ) THEN
+            ELSE IF( pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp ) THEN
                pt%xi = ricb_left + MOD(pt%xi, 1._wp )
             ENDIF
@@ -128,5 +126,5 @@
          pt => this%current_point
          ijne = INT( pt%yj + 0.5 )
-         IF( ijne .GT. mjg(nicbej) ) THEN
+         IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN
             !
             iine = INT( pt%xi + 0.5 )
@@ -170,5 +168,4 @@
       INTEGER                             ::   ibergs_rcvd_from_n, ibergs_rcvd_from_s
       INTEGER                             ::   i, ibergs_start, ibergs_end
-      INTEGER                             ::   iine, ijne
       INTEGER                             ::   ipe_N, ipe_S, ipe_W, ipe_E
       REAL(wp), DIMENSION(2)              ::   zewbergs, zwebergs, znsbergs, zsnbergs
@@ -234,6 +231,5 @@
          DO WHILE (ASSOCIATED(this))
             pt => this%current_point
-            iine = INT( pt%xi + 0.5 )
-            IF( ipe_E >= 0 .AND. iine > mig(nicbei) ) THEN
+            IF( ipe_E >= 0 .AND. pt%xi > REAL(mig(nicbei),wp) + 0.5_wp ) THEN
                tmpberg => this
                this => this%next
@@ -248,5 +244,5 @@
                CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e)
                CALL icb_utl_delete(first_berg, tmpberg)
-            ELSE IF( ipe_W >= 0 .AND. iine < mig(nicbdi) ) THEN
+            ELSE IF( ipe_W >= 0 .AND. pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp ) THEN
                tmpberg => this
                this => this%next
@@ -372,6 +368,5 @@
          DO WHILE (ASSOCIATED(this))
             pt => this%current_point
-            ijne = INT( pt%yj + 0.5 )
-            IF( ipe_N >= 0 .AND. ijne .GT. mjg(nicbej) ) THEN
+            IF( ipe_N >= 0 .AND. pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN
                tmpberg => this
                this => this%next
@@ -383,5 +378,5 @@
                CALL icb_pack_into_buffer( tmpberg, obuffer_n, ibergs_to_send_n)
                CALL icb_utl_delete(first_berg, tmpberg)
-            ELSE IF( ipe_S >= 0 .AND. ijne .LT. mjg(nicbdj) ) THEN
+            ELSE IF( ipe_S >= 0 .AND. pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp ) THEN
                tmpberg => this
                this => this%next
@@ -539,12 +534,10 @@
          DO WHILE (ASSOCIATED(this))
             pt => this%current_point
-            iine = INT( pt%xi + 0.5 )
-            ijne = INT( pt%yj + 0.5 )
-            IF( iine .LT. mig(nicbdi) .OR. &
-                iine .GT. mig(nicbei) .OR. &
-                ijne .LT. mjg(nicbdj) .OR. &
-                ijne .GT. mjg(nicbej)) THEN
+            IF( pt%xi < REAL(mig(nicbdi),wp) - 0.5_wp .OR. &
+                pt%xi > REAL(mig(nicbei),wp) + 0.5_wp .OR. &
+                pt%yj < REAL(mjg(nicbdj),wp) - 0.5_wp .OR. &
+                pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN
                i = i + 1
-               WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne
+               WRITE(numicb,*) 'berg lost in halo: ', this%number(:)
                WRITE(numicb,*) '                   ', nimpp, njmpp
                WRITE(numicb,*) '                   ', nicbdi, nicbei, nicbdj, nicbej
@@ -614,7 +607,6 @@
                   pt => this%current_point
                   iine = INT( pt%xi + 0.5 )
-                  ijne = INT( pt%yj + 0.5 )
                   iproc = nicbflddest(mi1(iine))
-                  IF( ijne .GT. mjg(nicbej) ) THEN
+                  IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN
                      IF( iproc == ifldproc ) THEN
                         !
@@ -696,5 +688,5 @@
                   ipts  = nicbfldpts (mi1(iine))
                   iproc = nicbflddest(mi1(iine))
-                  IF( ijne .GT. mjg(nicbej) ) THEN
+                  IF( pt%yj > REAL(mjg(nicbej),wp) + 0.5_wp ) THEN
                      IF( iproc == ifldproc ) THEN
                         !
Index: /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/LBC/lib_mpp.F90
===================================================================
--- /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/LBC/lib_mpp.F90	(revision 13180)
+++ /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/LBC/lib_mpp.F90	(revision 13181)
@@ -261,8 +261,14 @@
       !!
       INTEGER ::   iflag
-      !!----------------------------------------------------------------------
-      !
-#if defined key_mpp_mpi
-      CALL mpi_isend( pmess, kbytes, mpi_double_precision, kdest , ktyp, mpi_comm_oce, md_req, iflag )
+      INTEGER :: mpi_working_type
+      !!----------------------------------------------------------------------
+      !
+#if defined key_mpp_mpi
+      IF (wp == dp) THEN
+         mpi_working_type = mpi_double_precision
+      ELSE
+         mpi_working_type = mpi_real
+      END IF
+      CALL mpi_isend( pmess, kbytes, mpi_working_type, kdest , ktyp, mpi_comm_oce, md_req, iflag )
 #endif
       !
@@ -331,4 +337,5 @@
       INTEGER :: iflag
       INTEGER :: use_source
+      INTEGER :: mpi_working_type
       !!----------------------------------------------------------------------
       !
@@ -339,5 +346,10 @@
       IF( PRESENT(ksource) )   use_source = ksource
       !
-      CALL mpi_recv( pmess, kbytes, mpi_double_precision, use_source, ktyp, mpi_comm_oce, istatus, iflag )
+      IF (wp == dp) THEN
+         mpi_working_type = mpi_double_precision
+      ELSE
+         mpi_working_type = mpi_real
+      END IF
+      CALL mpi_recv( pmess, kbytes, mpi_working_type, use_source, ktyp, mpi_comm_oce, istatus, iflag )
 #endif
       !
Index: /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/grt_cir_dis.h90
===================================================================
--- /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/grt_cir_dis.h90	(revision 13180)
+++ /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/grt_cir_dis.h90	(revision 13181)
@@ -28,6 +28,12 @@
       REAL(KIND=wp) :: pc2   !  cos(lat2) * sin(lon2)
 
+      REAL(KIND=wp) :: cosdist ! cosine of great circle distance
+
+      ! Compute cosine of great circle distance, constraining it to be between
+      ! -1 and 1 (rounding errors can take it slightly outside this range
+      cosdist = MAX( MIN( pa1 * pa2 + pb1 * pb2 + pc1 * pc2, 1.0_wp), -1.0_wp )
+
       grt_cir_dis = &
-         &  ASIN( SQRT( 1.0 - ( pa1 * pa2 + pb1 * pb2 + pc1 * pc2 )**2 ) )
+         &  ASIN( SQRT( 1.0_wp - cosdist**2.0_wp ) )
       
    END FUNCTION grt_cir_dis
Index: /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obsinter_z1d.h90
===================================================================
--- /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obsinter_z1d.h90	(revision 13180)
+++ /NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obsinter_z1d.h90	(revision 13181)
@@ -62,29 +62,33 @@
          z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
          z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
-         IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp
+         
+         ! If kkco(jdep) is masked then set pobs(jdep) to the lowest value located above bathymetry
+         IF ( pobsmask(kkco(jdep)) == 0.0_wp ) THEN
+            pobs(jdep) = pobsk(kkco(jdep)-1)
+         ELSE
+            zsum = z1dm + z1dp
 
-         zsum = z1dm + z1dp
-         
-         IF ( k1dint == 0 ) THEN
+            IF ( k1dint == 0 ) THEN
 
-            !-----------------------------------------------------------------
-            !  Linear interpolation
-            !-----------------------------------------------------------------
-            pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
-               &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
+               !-----------------------------------------------------------------
+               !  Linear interpolation
+               !-----------------------------------------------------------------
+               pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
+                  &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
 
-         ELSEIF ( k1dint == 1 ) THEN
+            ELSEIF ( k1dint == 1 ) THEN
 
-            !-----------------------------------------------------------------
-            ! Cubic spline interpolation
-            !-----------------------------------------------------------------
-            zsum2 = zsum * zsum
-            pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &
-               &           + z1dp                             * pobsk (kkco(jdep)  ) &
-               &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &
-               &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &
-               &             ) / 6.0_wp                                              &
-               &          ) / zsum
+               !-----------------------------------------------------------------
+               ! Cubic spline interpolation
+               !-----------------------------------------------------------------
+               zsum2 = zsum * zsum
+               pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &
+                  &           + z1dp                             * pobsk (kkco(jdep)  ) &
+                  &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &
+                  &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &
+                  &             ) / 6.0_wp                                              &
+                  &          ) / zsum
 
+            ENDIF
          ENDIF
       END DO
