src/OCE/SBC/cpl_oasis3.F90

@@ -30,20 +30,35 @@ MODULE cpl_oasis3
    USE xios                         ! I/O server
 #endif
    USE par_oce                      ! ocean parameters
+
+   USE cpl_rnf_1d, ONLY: nn_cpl_river   ! Variables used in 1D river outflow 
+
    USE dom_oce                      ! ocean space and time domain
    USE in_out_manager               ! I/O manager
    USE lbclnk                       ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp
+#if defined key_agrif || ! defined key_mpi_off
+   USE MPI
+#endif
 
    IMPLICIT NONE
    PRIVATE
+#if ! defined key_oasis3
+   ! Dummy interface to oasis_get if not using oasis
+   INTERFACE oasis_get
+      MODULE PROCEDURE oasis_get_1d, oasis_get_2d
+   END INTERFACE
+#endif 
 
    PUBLIC   cpl_init
    PUBLIC   cpl_define
    PUBLIC   cpl_snd
    PUBLIC   cpl_rcv
+   PUBLIC   cpl_rcv_1d
    PUBLIC   cpl_freq
    PUBLIC   cpl_finalize
 
+   INTEGER, PARAMETER         ::   localRoot  = 0
    INTEGER, PUBLIC            ::   OASIS_Rcv  = 1    !: return code if received field
    INTEGER, PUBLIC            ::   OASIS_idle = 0    !: return code if nothing done by oasis
    INTEGER                    ::   ncomp_id          ! id returned by oasis_init_comp
@@ -67,9 +82,13 @@ MODULE cpl_oasis3
    INTEGER                    ::   nrcv         ! total number of fields received
    INTEGER                    ::   nsnd         ! total number of fields sent
    INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
-   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=62   ! Maximum number of coupling fields
+
+  INTEGER, PUBLIC, PARAMETER ::   nmaxfld=100  ! Maximum number of coupling fields
+
    INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields
+   LOGICAL, PARAMETER         ::   ltmp_wapatch = .FALSE.   ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
+   LOGICAL, PARAMETER         ::   ltmp_landproc = .TRUE.   ! patch to restrict coupled area to non halo cells
 
    TYPE, PUBLIC ::   FLD_CPL               !: Type for coupling field information
       LOGICAL               ::   laction   ! To be coupled or not
@@ -79,11 +98,16 @@ MODULE cpl_oasis3
       INTEGER, DIMENSION(nmaxcat,nmaxcpl) ::   nid   ! Id of the field (no more than 9 categories and 9 extrena models)
       INTEGER               ::   nct       ! Number of categories in field
       INTEGER               ::   ncplmodel ! Maximum number of models to/from which this variable may be sent/received
+      INTEGER               ::   dimensions ! Number of dimensions of coupling field 
    END TYPE FLD_CPL
 
    TYPE(FLD_CPL), DIMENSION(nmaxfld), PUBLIC ::   srcv, ssnd   !: Coupling fields
 
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   exfld   ! Temporary buffer for receiving
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   exfld_ext  ! Temporary buffer for receiving with wrap points
+
+   INTEGER :: ishape_ext(4)  ! shape of 2D arrays passed to PSMILe extended for wrap points in weights data
+   INTEGER :: ishape(4)      ! shape of 2D arrays passed to PSMILe 
 
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -103,6 +127,7 @@ CONTAINS
       !!--------------------------------------------------------------------
       CHARACTER(len = *), INTENT(in   ) ::   cd_modname   ! model name as set in namcouple file
       INTEGER           , INTENT(  out) ::   kl_comm      ! local communicator of the model
+      INTEGER                           ::   error
       !!--------------------------------------------------------------------
 
       ! WARNING: No write in numout in this routine
@@ -123,6 +148,7 @@ CONTAINS
       IF( nerror /= OASIS_Ok ) &
          CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' )
       !
+
    END SUBROUTINE cpl_init
 
 
@@ -138,13 +164,26 @@ CONTAINS
       INTEGER, INTENT(in) ::   krcv, ksnd     ! Number of received and sent coupling fields
       INTEGER, INTENT(in) ::   kcplmodel      ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
       !
-      INTEGER :: id_part
+      INTEGER :: id_part_0d     ! Partition for 0d fields
+      INTEGER :: id_part_rnf_1d ! Partition for 1d river outflow fields
+      INTEGER :: id_part_2d     ! Partition for 2d fields
+      INTEGER :: id_part_2d_ext ! Partition for 2d fields extended for old (pre vn4.2) style remapping weights!
+      INTEGER :: id_part_temp   ! Temperary partition used to choose either 0d or 1d partitions
       INTEGER :: paral(5)       ! OASIS3 box partition
-      INTEGER :: ishape(4)      ! shape of arrays passed to PSMILe
+
+      INTEGER :: paral_ext(5)       ! OASIS3 box partition extended
+      INTEGER :: ishape0d1d(2)  ! Shape of 0D or 1D arrays passed to PSMILe.
+      INTEGER :: var_nodims(2)  ! Number of coupling field dimensions.
+                                ! var_nodims(1) is redundant from OASIS3-MCT vn4.0 onwards
+                                ! but retained for backward compatibility.
+                                ! var_nodims(2) is the number of fields in a bundle
+                                ! or 1 for unbundled fields (bundles are not yet catered for
+                                ! in NEMO hence we default to 1).        
       INTEGER :: ji,jc,jm       ! local loop indicees
       LOGICAL :: llenddef       ! should we call xios_oasis_enddef and oasis_enddef?
       CHARACTER(LEN=64) :: zclname
       CHARACTER(LEN=2) :: cli2
+      INTEGER :: i_offset       ! Used in calculating offset for extended partition.
       !!--------------------------------------------------------------------
 
       IF(lwp) WRITE(numout,*)
@@ -173,10 +212,13 @@ CONTAINS
       ishape(2) = Ni_0
       ishape(3) = 1
       ishape(4) = Nj_0
-      !
+
+      ishape0d1d(1) = 0
+      ishape0d1d(2) = 0       !
+
       ! ... Allocate memory for data exchange
       !
-      ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror)        ! allocate only inner domain (without halos)
+      ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror)                 ! allocate full domain (without halos)
       IF( nerror > 0 ) THEN
          CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN
       ENDIF
@@ -198,7 +240,96 @@ CONTAINS
          WRITE(numout,*) ' multiexchg: Njs0, Nje0, njmpp =', Njs0, Nje0, njmpp
       ENDIF
 
-      CALL oasis_def_partition ( id_part, paral, nerror, Ni0glo*Nj0glo )   ! global number of points, excluding halos
+
+      ! We still set up the new vn4.2 style box partition for reference, though it doesn't actually get used,
+      ! we can easily swap back to it if we ever manage to successfully generate vn4.2 compatible weights, or introduce 
+      ! RTL controls to distinguish between onl and new style weights.  
+
+      CALL oasis_def_partition ( id_part_2d, paral, nerror, Ni0glo*Nj0glo )   ! global number of points, excluding halos
+
+      ! RSRH Set up 2D box partition for compatibility with existing weights files
+      ! so we don't have to generate and manage multiple sets of weights purely because of 
+      ! the changes to nemo 4.2+ code!
+
+      ! This is just a hack for global cyclic models for the time being
+      Ni0glo_ext = jpiglo
+      Nj0glo_ext = Nj0glo +1 ! We can't use jpjglo here because for some reason at 4.2 this is bigger
+                             ! than at 4.0.... e.g. for ORCA1 it is 333 when it should only be 332!
+
+      ! RSRH extended shapes for old style dimensioning. Allows backwards compatibility with existing weights files, 
+      ! which the new code DOES NOT, causing headaches not only for users but also for management of weights files. 
+      ishape_ext(:) = ishape(:)
+      IF (mig(1) == 1 .OR. mig(jpi)==jpiglo) THEN
+         ! Extra columns in PEs dealing with wrap points
+         ishape_ext(2) = ishape_ext(2) + 1
+      ENDIF
+
+      ! Workout any extra offset in the i dimension
+      IF (mig(1) == 1 ) THEN
+         i_offset = mig0(nn_hls)
+      ELSE
+         i_offset = mig(nn_hls)
+      ENDIF
+       
+      ALLOCATE(exfld_ext(ishape_ext(2), ishape_ext(4)), stat = nerror)        ! allocate full domain (with wrap pts)
+      IF( nerror > 0 ) THEN
+         CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld_ext')   ;   RETURN
+      ENDIF
+
+
+      ! Now we have the appropriate dimensions, we can set up the partition array for the old-style extended grid
+      paral_ext(1) = 2                                      ! box partitioning
+      paral_ext(2) = (Ni0glo_ext * mjg0(nn_hls)) + i_offset ! NEMO lower left corner global offset, with wrap pts
+      paral_ext(3) = Ni_0_ext                               ! local extent in i, including halos
+      paral_ext(4) = Nj_0_ext                               ! local extent in j, including halos
+      paral_ext(5) = Ni0glo_ext                             ! global extent in x, including halos
+
+      IF( sn_cfctl%l_oasout ) THEN
+         WRITE(numout,*) ' multiexchg: paral_ext (1:5)', paral_ext, jpiglo, jpjglo, Ni0glo_ext, Nj0glo_ext
+         WRITE(numout,*) ' multiexchg: Ni_0_ext, Nj_0_ext i_offset =', Ni_0_ext, Nj_0_ext, i_offset
+         WRITE(numout,*) ' multiexchg: Nis0_ext, Nie0_ext =', Nis0_ext, Nie0_ext
+         WRITE(numout,*) ' multiexchg: Njs0_ext, Nje0_ext =', Njs0_ext, Nje0_ext
+      ENDIF
+
+      ! Define our extended grid
+      CALL oasis_def_partition ( id_part_2d_ext, paral_ext, nerror, Ni0glo_ext*Nj0glo_ext ) 
+  
+      ! OK so now we should have a usable 2d partition for fields defined WITH redundant points. 
+
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      ! A special partition is needed for 0D fields
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     
+      paral(1) = 0                                       ! serial partitioning
+      paral(2) = 0   
+      IF ( mpprank == 0) THEN
+         paral(3) = 1                   ! Size of array to couple (scalar)
+      ELSE
+         paral(3) = 0                   ! Dummy size for PE's not involved
+      END IF
+      paral(4) = 0
+      paral(5) = 0
+             
+      CALL oasis_def_partition ( id_part_0d, paral, nerror, 1 )
+
+
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      ! Another special partition is needed for 1D river routing fields
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+     
+      paral(1) = 0                                       ! serial partitioning
+      paral(2) = 0   
+      IF ( mpprank == 0) THEN
+         paral(3) = nn_cpl_river                   ! Size of array to couple (vector)
+      ELSE
+         paral(3) = 0                   ! Dummy size for PE's not involved
+      END IF
+      paral(4) = 0
+      paral(5) = 0
+
+
+      CALL oasis_def_partition ( id_part_rnf_1d, paral, nerror, nn_cpl_river )
+
       !
       ! ... Announce send variables.
       !
@@ -232,14 +363,24 @@ CONTAINS
                   ENDIF
 #endif
                   IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
-                  CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
-                     &                OASIS_Out          , ishape , OASIS_REAL, nerror )
+
+                  flush(numout)
+
+
+                  IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out
+                     CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part_2d_ext   , (/ 2, 1 /),   &
+                     &                OASIS_Out           , ishape_ext , OASIS_REAL, nerror )
+
                   IF( nerror /= OASIS_Ok ) THEN
                      WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
                      CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
                   ENDIF
-                  IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "variable defined in the namcouple"
-                  IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "variable NOT defined in the namcouple"
+                  IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "put variable defined in the namcouple"
+                  IF( sn_cfctl%l_oasout .AND. ssnd(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "put variable NOT defined in the namcouple"
+
+
+
+
                END DO
             END DO
          ENDIF
@@ -276,15 +417,57 @@ CONTAINS
                      zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname)
                   ENDIF
 #endif
+
                   IF( sn_cfctl%l_oasout ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In
-                  CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   &
-                     &                OASIS_In           , ishape , OASIS_REAL, nerror )
+
+                  flush(numout)
+
+                  ! Define 0D (Greenland or Antarctic ice mass) or 1D (river outflow) coupling fields
+                  IF (srcv(ji)%dimensions <= 1) THEN
+                     var_nodims(1) = 1
+                     var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
+                     IF (mpprank == 0) THEN
+	
+                       IF (srcv(ji)%dimensions == 0) THEN
+
+                          ! If 0D then set temporary variables to 0D components
+                          id_part_temp = id_part_0d
+                          ishape0d1d(2) = 1
+                       ELSE
+
+                          ! If 1D then set temporary variables to river outflow components
+                          id_part_temp = id_part_rnf_1d
+                          ishape0d1d(2)= nn_cpl_river
+
+                       END IF
+
+                       CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_temp   , var_nodims,   &
+                                        OASIS_In           , ishape0d1d(1:2) , OASIS_REAL, nerror )
+                    ELSE
+                       ! Dummy call to keep OASIS3-MCT happy.
+                       CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_0d   , var_nodims,   &
+                                        OASIS_In           , ishape0d1d(1:2) , OASIS_REAL, nerror )
+                    END IF
+                  ELSE
+                    ! It's a "normal" 2D (or pseudo 3D) coupling field.
+                    ! ... Set the field dimension and bundle count
+                    var_nodims(1) = 2
+                    var_nodims(2) = 1 ! Modify this value to cater for bundled fields.
+
+                    CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part_2d_ext   , var_nodims,   &
+                                        OASIS_In           , ishape_ext , OASIS_REAL, nerror )
+
+                  ENDIF 
+
+
+
                   IF( nerror /= OASIS_Ok ) THEN
                      WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname)
                      CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' )
                   ENDIF
-                  IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "variable defined in the namcouple"
-                  IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "variable NOT defined in the namcouple"
+                  IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) /= -1 ) WRITE(numout,*) "get variable defined in the namcouple"
+                  IF( sn_cfctl%l_oasout .AND. srcv(ji)%nid(jc,jm) == -1 ) WRITE(numout,*) "get variable NOT defined in the namcouple"
+
 
                END DO
             END DO
@@ -330,13 +513,17 @@ CONTAINS
       INTEGER                                   ::   jc,jm     ! local loop index
       !!--------------------------------------------------------------------
       !
+
       ! snd data to OASIS3
       !
       DO jc = 1, ssnd(kid)%nct
          DO jm = 1, ssnd(kid)%ncplmodel
 
             IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN   ! exclude halos from data sent to oasis
-               CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0:Nie0, Njs0:Nje0,jc), kinfo )
+
+               ! The field is "put" directly, using appropriate start/end indexing - i.e. we don't
+               ! copy it to an intermediate buffer. 
+               CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc), kinfo )
 
                IF ( sn_cfctl%l_oasout ) THEN
                   IF ( kinfo == OASIS_Sent     .OR. kinfo == OASIS_ToRest .OR.   &
@@ -346,10 +533,11 @@ CONTAINS
                      WRITE(numout,*) 'oasis_put: ivarid ', ssnd(kid)%nid(jc,jm)
                      WRITE(numout,*) 'oasis_put:  kstep ', kstep
                      WRITE(numout,*) 'oasis_put:   info ', kinfo
-                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
-                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
-                     WRITE(numout,*) '     -     Sum value is ',    SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc))
+                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
+                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
+                     WRITE(numout,*) '     -     Sum value is ',    SUM(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
                      WRITE(numout,*) '****************'
+                     CALL FLUSH(numout)
                   ENDIF
                ENDIF
 
@@ -357,6 +545,7 @@ CONTAINS
 
          ENDDO
       ENDDO
+
       !
     END SUBROUTINE cpl_snd
 
@@ -375,13 +564,16 @@ CONTAINS
       INTEGER                   , INTENT(  out) ::   kinfo     ! OASIS3 info argument
       !!
       INTEGER                                   ::   jc,jm     ! local loop index
-      LOGICAL                                   ::   llaction, ll_1st
+      LOGICAL                                   ::   llaction, ll_1st, lrcv
+
       !!--------------------------------------------------------------------
       !
       ! receive local data from OASIS3 on every process
       !
       kinfo = OASIS_idle
+      lrcv=.FALSE.
       !
+
       DO jc = 1, srcv(kid)%nct
          ll_1st = .TRUE.
 
@@ -389,7 +581,7 @@ CONTAINS
 
             IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
 
-               CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld, kinfo )
+               CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, exfld_ext, kinfo )
 
                llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   &
                   &        kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
@@ -398,14 +590,18 @@ CONTAINS
                   &  WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
 
                IF( llaction ) THEN   ! data received from oasis do not include halos
+                                     ! but DO still cater for wrap columns when using pre vn4.2 compatible remapping weights. 
 
                   kinfo = OASIS_Rcv
+                  lrcv=.TRUE. 
                   IF( ll_1st ) THEN
-                     pdata(Nis0:Nie0,Njs0:Nje0,jc) =   exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
+                     pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) =   exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
+
                      ll_1st = .FALSE.
                   ELSE
-                     pdata(Nis0:Nie0,Njs0:Nje0,jc) = pdata(Nis0:Nie0,Njs0:Nje0,jc)   &
-                        &                                + exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm)
+
+                     pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc) = pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc)   &
+                        &                                + exfld_ext(:,:) * pmask(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jm)
                   ENDIF
 
                   IF ( sn_cfctl%l_oasout ) THEN
@@ -414,10 +610,11 @@ CONTAINS
                      WRITE(numout,*) 'oasis_get: ivarid '  , srcv(kid)%nid(jc,jm)
                      WRITE(numout,*) 'oasis_get:   kstep', kstep
                      WRITE(numout,*) 'oasis_get:   info ', kinfo
-                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
-                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc))
-                     WRITE(numout,*) '     -     Sum value is ',    SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc))
+                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
+                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
+                     WRITE(numout,*) '     -     Sum value is ',    SUM(pdata(Nis0_ext:Nie0_ext,Njs0_ext:Nje0_ext,jc))
                      WRITE(numout,*) '****************'
+                     CALL FLUSH(numout)
                   ENDIF
 
                ENDIF
@@ -426,15 +623,130 @@ CONTAINS
 
          ENDDO
 
-         !--- we must call lbc_lnk to fill the halos that where not received.
-         IF( .NOT. ll_1st ) THEN
-            CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )
-         ENDIF
-
       ENDDO
+
+      ! RSRH I've changed this since:
+      ! 1) it seems multi cat fields may not be properly updated in halos when called on a per 
+      !    category basis(?) 
+      ! 2) it's more efficient to have a single call (and simpler to understand) to update ALL 
+      !    categories at the same time!
+      !--- Call lbc_lnk to populate halos of received fields.
+      IF (lrcv) then
+         CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,:), srcv(kid)%clgrid, srcv(kid)%nsgn )
+      endif
       !
    END SUBROUTINE cpl_rcv
 
+  SUBROUTINE cpl_rcv_1d( kid, kstep, pdata, nitems, kinfo )
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE cpl_rcv_1d  ***
+      !!
+      !! ** Purpose : - A special version of cpl_rcv to deal exclusively with
+      !! receipt of 0D or 1D fields.
+      !! The fields are recieved into a 1D array buffer which is simply a
+      !! dynamically sized sized array (which may be of size 1)
+      !! of 0 dimensional fields. This allows us to pass miltiple 0D
+      !! fields via a single put/get operation. 
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in   ) ::   nitems      ! Number of 0D items to recieve
+                                               ! during this get operation. i.e.
+                                               ! The size of the 1D array in which
+                                               ! 0D items are passed.   
+      INTEGER , INTENT(in   ) ::   kid         ! ID index of the incoming
+                                               ! data. 
+      INTEGER , INTENT(in   ) ::   kstep       ! ocean time-step in seconds
+      REAL(wp), INTENT(inout) ::   pdata(1:nitems) ! The original value(s), 
+                                                   ! unchanged if nothing is
+                                                   ! received
+      INTEGER , INTENT(  out) ::   kinfo       ! OASIS3 info argument
+      !!
+      REAL(wp) ::   recvfld(1:nitems)          ! Local receive field buffer
+      INTEGER  ::   jc,jm     ! local loop index
+      INTEGER  ::   ierr
+      LOGICAL  ::   llaction
+      INTEGER  ::   MPI_WORKING_PRECISION
+      INTEGER  ::   number_to_print
+      !!--------------------------------------------------------------------
+      !
+      ! receive local data from OASIS3 on every process
+      !
+      kinfo = OASIS_idle
+      !
+      ! 0D and 1D fields won't have categories or any other form of "pseudo level"
+      ! so we only cater for a single set of values and thus don't bother
+      ! with a loop over the jc index
+      jc = 1
+
+      DO jm = 1, srcv(kid)%ncplmodel
+
+         IF( srcv(kid)%nid(jc,jm) /= -1 ) THEN
+
+            IF ( ( srcv(kid)%dimensions <= 1) .AND. (mpprank == 0) ) THEN
+                ! Since there is no concept of data decomposition for zero
+                ! dimension fields, they must only be exchanged through the master PE,
+                ! unlike "normal" 2D field cases where every PE is involved.
+ 
+                CALL oasis_get ( srcv(kid)%nid(jc,jm), kstep, recvfld, kinfo )   
+                
+                llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   &
+                            kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut
+                
+                IF ( sn_cfctl%l_oasout ) WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , &
+                                      llaction, kinfo, kstep, srcv(kid)%nid(jc,jm)
+                
+                IF ( llaction ) THEN
+	                 
+                   kinfo = OASIS_Rcv
+	           pdata(1:nitems) = recvfld(1:nitems)
+            
+                   IF( sn_cfctl%l_oasout ) THEN     
+	                     number_to_print = 10
+	                     IF ( nitems < number_to_print ) number_to_print = nitems
+	                     WRITE(numout,*) '****************'
+	                     WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname
+	                     WRITE(numout,*) 'oasis_get: ivarid '  , srcv(kid)%nid(jc,jm)
+	                     WRITE(numout,*) 'oasis_get:   kstep', kstep
+	                     WRITE(numout,*) 'oasis_get:   info ', kinfo
+	                     WRITE(numout,*) '     - Minimum Value is ', MINVAL(pdata(:))
+	                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(:))
+	                     WRITE(numout,*) '     - Start of data is ', pdata(1:number_to_print)
+	                     WRITE(numout,*) '****************'
+                             CALL FLUSH(numout)
+                  ENDIF
+	                 
+               ENDIF
+            ENDIF   
+          ENDIF
+	           
+       ENDDO
+
+#if defined key_mpi_off
+
+       CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Unable to use mpi_bcast with key_mpi_off in your list of NEMO keys." )
+
+#else
+       ! Set the precision that we want to broadcast using MPI_BCAST
+       SELECT CASE( wp )
+       CASE( sp )
+         MPI_WORKING_PRECISION = MPI_REAL                ! Single precision
+       CASE( dp )
+         MPI_WORKING_PRECISION = MPI_DOUBLE_PRECISION    ! Double precision
+       CASE default
+         CALL oasis_abort( ncomp_id, "cpl_rcv_1d", "Could not find precision for coupling 0d or 1d field" )
+       END SELECT
+
+       ! We have to broadcast (potentially) received values from PE 0 to all
+       ! the others. If no new data has been received we're just
+       ! broadcasting the existing values but there's no more efficient way
+       ! to deal with that w/o NEMO adopting a UM-style test mechanism
+       ! to determine active put/get timesteps.
+       CALL mpi_bcast( pdata, nitems, MPI_WORKING_PRECISION, localRoot, mpi_comm_oce, ierr )
+
+#endif
+	
+        !
+END SUBROUTINE cpl_rcv_1d
+
 
    INTEGER FUNCTION cpl_freq( cdfieldname )
       !!---------------------------------------------------------------------
@@ -500,6 +812,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       !
       DEALLOCATE( exfld )
+      DEALLOCATE( exfld_ext )
       IF(nstop == 0) THEN
          CALL oasis_terminate( nerror )
       ELSE
@@ -543,7 +856,7 @@ CONTAINS
 
    SUBROUTINE oasis_def_var(k1,cd1,k2,k3,k4,k5,k6,k7)
       CHARACTER(*), INTENT(in   ) ::  cd1
-      INTEGER     , INTENT(in   ) ::  k2,k3(2),k4,k5(2,2),k6
+      INTEGER     , INTENT(in   ) ::  k2,k3(2),k4,k5(*),k6
       INTEGER     , INTENT(  out) ::  k1,k7
       k1 = -1 ; k7 = -1
       WRITE(numout,*) 'oasis_def_var: Error you sould not be there...', cd1
@@ -563,13 +876,21 @@ CONTAINS
       WRITE(numout,*) 'oasis_put: Error you sould not be there...'
    END SUBROUTINE oasis_put
 
-   SUBROUTINE oasis_get(k1,k2,p1,k3)
+   SUBROUTINE oasis_get_2d(k1,k2,p1,k3)
       REAL(wp), DIMENSION(:,:), INTENT(  out) ::  p1
       INTEGER                 , INTENT(in   ) ::  k1,k2
       INTEGER                 , INTENT(  out) ::  k3
       p1(1,1) = -1. ; k3 = -1
-      WRITE(numout,*) 'oasis_get: Error you sould not be there...'
-   END SUBROUTINE oasis_get
+      WRITE(numout,*) 'oasis_get_2d: Error you sould not be there...'
+   END SUBROUTINE oasis_get_2d
+
+   SUBROUTINE oasis_get_1d(k1,k2,p1,k3)
+      REAL(wp), DIMENSION(:)  , INTENT(  out) ::  p1
+      INTEGER                 , INTENT(in   ) ::  k1,k2
+      INTEGER                 , INTENT(  out) ::  k3
+      p1(1) = -1. ; k3 = -1
+      WRITE(numout,*) 'oasis_get_1d: Error you sould not be there...'
+   END SUBROUTINE oasis_get_1d
 
    SUBROUTINE oasis_get_freqs(k1,k5,k2,k3,k4)
       INTEGER              , INTENT(in   ) ::  k1,k2

src/OCE/SBC/cpl_rnf_1d.F90

+MODULE cpl_rnf_1d
+   !!======================================================================
+   !!                       ***  MODULE  cpl_rnf_1d  ***
+   !! Ocean forcing:  River runoff passed from the atmosphere using
+   !!                 1D array. One value per river.
+   !!=====================================================================
+   !! History : ?.?  ! 2018-01 (D. Copsey) Initial setup
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   cpl_rnf_1d_init : runoffs initialisation
+   !!---------------------------------------------------------------------- 
+
+#if defined key_oasis3
+   USE mod_oasis                    ! OASIS3-MCT module
+#endif
+   USE timing          ! Timing
+   USE in_out_manager  ! I/O units
+   USE lib_mpp         ! MPP library
+   USE iom
+   USE dom_oce         ! Domain sizes (for grid box area e1e2t)
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE lib_fortran,    ONLY: DDPDD
+   
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   cpl_rnf_1d_init     ! routine called in nemo_init
+   PUBLIC   cpl_rnf_1d_to_2d      ! routine called in sbccpl.F90
+   
+   TYPE, PUBLIC ::   RIVERS_DATA     !: Storage for river outflow data
+      INTEGER, ALLOCATABLE, DIMENSION(:,:)    ::   river_number       !: River outflow number
+      REAL(wp), ALLOCATABLE, DIMENSION(:)     ::   river_area         ! 1D array listing areas of each river outflow (m2)
+      COMPLEX(wp), ALLOCATABLE, DIMENSION(:)  ::   river_area_c       ! Comlex version of river_area for use in bit reproducible sums (m2)
+   END TYPE RIVERS_DATA
+   
+   ! DB: see https://forge.ipsl.jussieu.fr/nemo/ticket/1865
+#if defined key_agrif
+   TYPE(RIVERS_DATA), PUBLIC         :: rivers  !: River data
+#else
+   TYPE(RIVERS_DATA), PUBLIC, TARGET :: rivers  !: River data
+#endif
+
+   INTEGER, PUBLIC            :: nn_cpl_river   ! Maximum number of rivers being passed through the coupler
+   INTEGER, PUBLIC            :: runoff_id      ! OASIS coupling id used in oasis_get command
+   LOGICAL                    :: ln_print_river_info  ! Diagnostic prints of river coupling information
+   
+CONTAINS
+
+   SUBROUTINE cpl_rnf_1d_init
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE cpl_rnf_1d_init  ***
+      !!                     
+      !! ** Purpose : - Read in file for river outflow numbers.
+      !!                Calculate 2D area of river outflow points.
+      !!                Called from nemo_init (nemogcm.F90).
+      !!
+      !!----------------------------------------------------------------------
+      !! namelist variables
+      !!-------------------
+      CHARACTER(len=200)                        ::   file_riv_number             !: Filename for river numbers
+      INTEGER                                   ::   ios                 ! Local integer output status for namelist read
+      INTEGER                                   ::   inum
+      INTEGER                                   ::   ii, jj, rr          !: Loop indices
+      INTEGER                                   ::   max_river
+      REAL(wp), DIMENSION(jpi,jpj)              ::   river_number        ! 2D array containing the river outflow numbers
+      
+      NAMELIST/nam_cpl_rnf_1d/file_riv_number, nn_cpl_river, ln_print_river_info
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing ) CALL timing_start('cpl_rnf_1d_init')
+      
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'cpl_rnf_1d_init : initialization of river runoff coupling'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      IF(lwp) CALL flush(numout)
+      
+!!      REWIND(numnam_cfg)
+      
+      ! Read the namelist
+      READ  ( numnam_ref, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in reference namelist' )
+      READ  ( numnam_cfg, nam_cpl_rnf_1d, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_cpl_rnf_1d in configuration namelist' )
+      IF(lwm) WRITE ( numond, nam_cpl_rnf_1d )
+
+      !                                               ! Parameter control and print
+      IF(lwp) WRITE(numout,*) '  '
+      IF(lwp) WRITE(numout,*) '          Namelist nam_cpl_rnf_1d : Coupled runoff using 1D array'
+      IF(lwp) WRITE(numout,*) '             Input file that contains river numbers = ',file_riv_number
+      IF(lwp) WRITE(numout,*) '             Maximum number of rivers to couple = ',nn_cpl_river
+      IF(lwp) WRITE(numout,*) '             Print river information = ',ln_print_river_info
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) CALL flush(numout)
+
+      ! Assign space for river numbers
+      ALLOCATE( rivers%river_number( jpi, jpj ) )
+      
+      ! Read the river numbers from netcdf file
+      CALL iom_open (file_riv_number , inum )
+      CALL iom_get  ( inum, jpdom_global, 'river_number', river_number, cd_type = 'T', psgn = 1._wp)
+      CALL iom_close( inum )
+      
+      ! Convert from a real array to an integer array
+      max_river=0
+      DO ii = 1, jpi
+        DO jj = 1, jpj
+          rivers%river_number(ii,jj) = INT(river_number(ii,jj))
+          IF ( rivers%river_number(ii,jj) > max_river ) THEN
+            max_river = rivers%river_number(ii,jj)
+          END IF
+        END DO
+      END DO
+      
+      ! Print out the largest river number
+      IF ( ln_print_river_info .AND. lwp) THEN
+         WRITE(numout,*) 'Maximum river number in input file = ',max_river
+         CALL flush(numout)
+      END IF
+      
+      ! Get the area of each river outflow
+      ALLOCATE( rivers%river_area( nn_cpl_river ) )
+      ALLOCATE( rivers%river_area_c( nn_cpl_river ) )
+      rivers%river_area_c(:) = CMPLX( 0.e0, 0.e0, wp )
+      DO ii = Nis0, Nie0     
+        DO jj = Njs0, Nje0
+          IF ( tmask_i(ii,jj) > 0.5 ) THEN  ! This makes sure we are not at a duplicated point (at north fold or east-west cyclic point)
+            IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN
+              ! Add the area of each grid box (e1e2t) into river_area_c using DDPDD which should maintain bit reproducibility (needs to be checked)
+              CALL DDPDD( CMPLX( e1e2t(ii,jj), 0.e0, wp ), rivers%river_area_c(rivers%river_number(ii,jj) ) )
+            END IF
+          END IF
+        END DO
+      END DO
+      
+      ! Use mpp_sum to add together river areas on other processors
+      CALL mpp_sum( 'cpl_rnf_1d', rivers%river_area_c )
+      
+      ! Convert from complex number to real
+      rivers%river_area(:) = REAL(rivers%river_area_c(:),wp)
+      
+      IF ( ln_print_river_info .AND. lwp) THEN
+        WRITE(numout,*) 'Area of rivers 1 to 10 are ',rivers%river_area(1:10)
+        WRITE(numout,*) 'Maximum river area = ',MAXVAL(rivers%river_area)
+        WRITE(numout,*) 'Minimum river area = ',MINVAL(rivers%river_area)
+        WRITE(numout,*) 'Sum of river areas = ',SUM(rivers%river_area)
+        CALL flush(numout)
+      END IF
+
+      IF ( MINVAL(rivers%river_area) <= 0 ) THEN
+         WRITE(numout,*) 'ERROR: There is at least one river with a river outflow area of zero. Please check your file_riv_number file'
+         WRITE(numout,*) 'that all the allocated river numbers are at ocean points as defined by the bathymetry file with no river'
+         WRITE(numout,*) 'numbers within the north fold or wraparound points.'
+         DO rr = 1,nn_cpl_river
+           IF ( rivers%river_area(rr) <= 0 ) THEN
+             WRITE(numout,*) 'The first river with this problem is river number ',rr
+             CALL ctl_stop ( 'STOP', 'ERROR: There is at least one river with a river outflow area of zero. See stdout.')
+           END IF
+         END DO
+      END IF
+      
+   END SUBROUTINE cpl_rnf_1d_init
+   
+   SUBROUTINE cpl_rnf_1d_to_2d( runoff_1d )
+   
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE cpl_rnf_1d_to_2d  ***
+      !!                     
+      !! ** Purpose : - Convert river outflow from 1D array (passed from the
+      !!                atmosphere) to the 2D NEMO runoff field.
+      !!                Called from sbc_cpl_ice_flx (sbccpl.F90).
+      !!
+      !!----------------------------------------------------------------------
+      
+      REAL                   , INTENT(in   ) ::   runoff_1d(nn_cpl_river)    ! River runoff. One value per river.
+      
+      INTEGER  ::   ii, jj, rr                 ! Loop indices
+      LOGICAL  ::   found_nan
+            
+      ! Convert the 1D total runoff per river to 2D runoff flux by
+      ! dividing by the area of each runoff zone.
+      DO ii = 1, jpi
+        DO jj = 1, jpj
+          IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river .AND. tmask_i(ii,jj) > 0.5 ) THEN
+            rnf(ii,jj) = runoff_1d(rivers%river_number(ii,jj)) / rivers%river_area(rivers%river_number(ii,jj))
+          ELSE
+            rnf(ii,jj) = 0.0
+          END IF
+            
+        END DO
+      END DO
+
+      IF ( ln_print_river_info ) THEN
+        WRITE(numout,*) 'SUM of river runoff on 1D array = ',SUM(runoff_1d)
+        WRITE(numout,*) 'SUM of river runoff on 2D array = ',SUM(rnf)
+        found_nan = .false.
+        DO rr = 1, nn_cpl_river
+          IF (rr <= 10) THEN
+            WRITE(numout,*) 'River number ',rr,' has total runoff=',runoff_1d(rr),' and area=',rivers%river_area(rr),' making runoff flux=',runoff_1d(rr)/rivers%river_area(rr)
+          END IF
+
+          IF (runoff_1d(rr) /= runoff_1d(rr)) THEN
+            IF (.NOT. found_nan) THEN
+              WRITE(numout,*) 'WARNING: NAN found at river number ',rr
+              found_nan = .true.
+            END IF
+          END IF
+
+        END DO
+      END IF     
+
+
+   END SUBROUTINE cpl_rnf_1d_to_2d
+
+END MODULE cpl_rnf_1d

src/OCE/SBC/fldread.F90

@@ -67,6 +67,7 @@ MODULE fldread
       REAL(wp)                        ::   freqh        ! frequency of each flux file
       CHARACTER(len = 34)             ::   clvar        ! generic name of the variable in the NetCDF flux file
       LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F)
+      REAL(wp)                        ::   rec_shft     ! record shift (from -0.5 to +0.5) when ln_tint=T
       LOGICAL                         ::   ln_clim      ! climatology or not (T/F)
       CHARACTER(len = 8)              ::   clftyp       ! type of data file 'daily', 'monthly' or yearly'
       CHARACTER(len = 1)              ::   cltype       ! nature of grid-points: T, U, V...
@@ -75,7 +76,7 @@ MODULE fldread
       INTEGER , DIMENSION(2,2)        ::   nrec         ! before/after record (1: index, 2: second since Jan. 1st 00h of yr nit000)
       INTEGER                         ::   nbb          ! index of before values
       INTEGER                         ::   naa          ! index of after  values
-      INTEGER , ALLOCATABLE, DIMENSION(:) ::   nrecsec   ! 
+      INTEGER , ALLOCATABLE, DIMENSION(:) ::   nrecsec  ! records time in seconds
       REAL(wp), POINTER, DIMENSION(:,:,:  ) ::   fnow   ! input fields interpolated to now time step
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   fdta   ! 2 consecutive record of input fields
       CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
@@ -129,6 +130,7 @@ MODULE fldread
 
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "single_precision_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -210,7 +212,7 @@ CONTAINS
             ibb = sd(jf)%nbb   ;   iaa = sd(jf)%naa
             !
             IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
-               IF(lwp .AND. ( kt - nit000 <= 20 .OR. nitend - kt <= 20 ) ) THEN 
+               IF(lwp .AND. ( kt - nit000 <= 30 .OR. nitend - kt <= 30 ) ) THEN 
                   clfmt = "('   fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
                   WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &            
@@ -222,7 +224,7 @@ CONTAINS
                ztintb =  1. - ztinta
                sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,ibb) + ztinta * sd(jf)%fdta(:,:,:,iaa)
             ELSE   ! nothing to do...
-               IF(lwp .AND. ( kt - nit000 <= 20 .OR. nitend - kt <= 20 ) ) THEN
+               IF(lwp .AND. ( kt - nit000 <= 30 .OR. nitend - kt <= 30 ) ) THEN
                   clfmt = "('   fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
                   WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
@@ -755,7 +757,7 @@ CONTAINS
       INTEGER  :: idaysec               ! number of seconds in 1 day = NINT(rday)
       INTEGER  :: iyr, imt, idy, isecwk
       INTEGER  :: indexyr, indexmt
-      INTEGER  :: ireclast
+      INTEGER  :: ireclast, irecshft
       INTEGER  :: ishift, istart
       INTEGER, DIMENSION(2)  :: isave
       REAL(wp) :: zfreqs
@@ -791,6 +793,7 @@ CONTAINS
       ! previous file parameters
       IF( llprev ) THEN
          IF( sdjf%clftyp(1:4) == 'week'    ) THEN   ! find the day of the beginning of previous week
+            isecwk = ksec_week( sdjf%clftyp(6:8) )  ! seconds between the beginning of the week and half of current time step
             isecwk = isecwk + 7 * idaysec           ! seconds between the beginning of previous week and half of the time step
             llprevmt = isecwk > nsec_month          ! longer time since beginning of the previous week than the current month
             llprevyr = llprevmt .AND. nmonth == 1
@@ -811,8 +814,9 @@ CONTAINS
       ! next file parameters
       IF( llnext ) THEN
          IF( sdjf%clftyp(1:4) == 'week'    ) THEN   ! find the day of the beginning of next week
+            isecwk = ksec_week( sdjf%clftyp(6:8) )  ! seconds between the beginning of the week and half of current time step
             isecwk = 7 * idaysec - isecwk           ! seconds between half of the time step and the beginning of next week
-            llnextmt = isecwk > ( nmonth_len(nmonth)*idaysec - nsec_month )   ! larger than the seconds to the end of the month
+            llnextmt = isecwk >= ( nmonth_len(nmonth)*idaysec - nsec_month )   ! larger than the seconds to the end of the month
             llnextyr = llnextmt .AND. nmonth == 12
             iyr = nyear  + COUNT((/llnextyr/))
             imt = nmonth + COUNT((/llnextmt/)) - 12 * COUNT((/llnextyr/))
@@ -885,8 +889,9 @@ CONTAINS
       !
       IF( sdjf%ln_tint ) THEN   ! record time defined in the middle of the record, computed using an implementation
                                 ! of the rounded average that is valid over the full integer range
+         irecshft = NINT( sdjf%rec_shft * (sdjf%nrecsec(1) - sdjf%nrecsec(0)) )
          sdjf%nrecsec(1:sdjf%nreclast) = sdjf%nrecsec(0:sdjf%nreclast-1) / 2 + sdjf%nrecsec(1:sdjf%nreclast) / 2 + &
-            & MAX( MOD( sdjf%nrecsec(0:sdjf%nreclast-1), 2 ), MOD( sdjf%nrecsec(1:sdjf%nreclast), 2 ) )
+            & MAX( MOD( sdjf%nrecsec(0:sdjf%nreclast-1), 2 ), MOD( sdjf%nrecsec(1:sdjf%nreclast), 2 ) ) + irecshft
       END IF
       !
       sdjf%clname = fld_filename( sdjf, idy, imt, iyr )
@@ -966,6 +971,16 @@ CONTAINS
          sdf(jf)%freqh      = sdf_n(jf)%freqh
          sdf(jf)%clvar      = sdf_n(jf)%clvar
          sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint
+         sdf(jf)%rec_shft   = 0._wp
+         IF( sdf(jf)%ln_tint ) THEN
+            IF( sdf(jf)%clvar(1:1) == '+' ) THEN
+               sdf(jf)%rec_shft = +0.5_wp
+               sdf(jf)%clvar = sdf(jf)%clvar(2:)
+            ELSE IF( sdf(jf)%clvar(1:1) == '-' ) THEN
+               sdf(jf)%rec_shft = -0.5_wp
+               sdf(jf)%clvar = sdf(jf)%clvar(2:)
+            ENDIF
+         ENDIF
          sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim
          sdf(jf)%clftyp     = sdf_n(jf)%clftyp
          sdf(jf)%cltype     = 'T'   ! by default don't do any call to lbc_lnk in iom_get
@@ -1100,7 +1115,7 @@ CONTAINS
       CHARACTER (len=5) ::   clname   !
       INTEGER , DIMENSION(4) ::   ddims
       INTEGER                ::   isrc
-      REAL(wp), DIMENSION(jpi,jpj) ::   data_tmp
+      REAL(dp), DIMENSION(jpi,jpj) ::   data_tmp
       !!----------------------------------------------------------------------
       !
       IF( nxt_wgt > tot_wgts ) THEN

src/OCE/SBC/geo2ocean.F90

@@ -25,6 +25,7 @@ MODULE geo2ocean
    IMPLICIT NONE
    PRIVATE
 
+   PUBLIC   repcmo    ! called in sbccpl
    PUBLIC   rot_rep   ! called in sbccpl, fldread, and cyclone
    PUBLIC   geo2oce   ! called in sbccpl
    PUBLIC   oce2geo   ! called in sbccpl
@@ -50,6 +51,48 @@ MODULE geo2ocean
    !!----------------------------------------------------------------------
 CONTAINS
 
+
+   SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
+                       px2 , py2 , kchoix  )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE repcmo  ***
+      !!
+      !! ** Purpose :   Change vector componantes from a geographic grid to a
+      !!      stretched coordinates grid.
+      !!
+      !! ** Method  :   Initialization of arrays at the first call.
+      !!
+      !! ** Action  : - px2 : first  componante (defined at u point)
+      !!              - py2 : second componante (defined at v point)
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxu1, pyu1   ! geographic vector componantes at u-point
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pxv1, pyv1   ! geographic vector componantes at v-point
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   px2          ! i-componante (defined at u-point)
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( IN ) ::   &
+         kchoix   ! type of transformation
+                  ! = 1 change from geographic to model grid.
+                  ! =-1 change from model to geographic grid
+      !!----------------------------------------------------------------------
+ 
+      SELECT CASE (kchoix)
+      CASE ( 1)
+        ! Change from geographic to stretched coordinate
+        ! ----------------------------------------------
+     
+        CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
+        CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
+      CASE (-1)
+        ! Change from stretched to geographic coordinate
+        ! ----------------------------------------------
+     
+        CALL rot_rep( pxu1, pyu1, 'U', 'ij->e',px2 )
+        CALL rot_rep( pxv1, pyv1, 'V', 'ij->n',py2 )
+      END SELECT
+     
+   END SUBROUTINE repcmo	
+
    SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE rot_rep  ***

src/OCE/SBC/sbc_ice.F90

@@ -94,6 +94,7 @@ MODULE sbc_ice
    ! already defined in ice.F90 for SI3
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  h_i, h_s
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_i_last_couple   !: Sea ice fraction on categories at the last coupling point
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K]
 #endif

src/OCE/SBC/sbc_oce.F90

@@ -126,11 +126,13 @@ MODULE sbc_oce
    !!
    !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tprecip           !: total precipitation                          [Kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
+   REAL(dp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra         !: cloud cover (fraction of cloud in a gridcell) [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   greenland_icesheet_mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   antarctica_icesheet_mask 
 
    !!---------------------------------------------------------------------
    !! ABL Vertical Domain size
@@ -153,6 +155,33 @@ MODULE sbc_oce
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
 
+   !!----------------------------------------------------------------------
+   !!  Surface scalars of total ice sheet mass for Greenland and Antarctica,
+   !! passed from atmosphere to be converted to dvol and hence a freshwater
+   !! flux  by using old values. New values are saved in the dump, to become
+   !! old values next coupling timestep. Freshwater fluxes split between
+   !! sub iceshelf melting and iceberg calving, scalled to flux per second
+   !!----------------------------------------------------------------------
+   
+   REAL(wp), PUBLIC  :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed
+   REAL(wp), PUBLIC  :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed
+
+   ! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to
+   ! avoid circular dependencies.
+   INTEGER, PUBLIC     ::   nn_coupled_iceshelf_fluxes     ! =0 : total freshwater input from iceberg calving and ice shelf basal melting
+                                                           ! taken from climatologies used (no action in coupling routines).
+                                                           ! =1 :  use rate of change of mass of Greenland and Antarctic icesheets to set the
+                                                           ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+                                                           ! =2 :  specify constant freshwater inputs in this namelist to set the combined
+                                                           ! magnitude of iceberg calving and iceshelf melting freshwater fluxes.
+   LOGICAL, PUBLIC     ::   ln_iceshelf_init_atmos         ! If true force ocean to initialise iceshelf masses from atmospheric values rather
+                                                           ! than values in ocean restart (applicable if nn_coupled_iceshelf_fluxes=1).
+   REAL(wp), PUBLIC    ::   rn_greenland_total_fw_flux    ! Constant total rate of freshwater input (kg/s) for Greenland (if nn_coupled_iceshelf_fluxes=2)
+   REAL(wp), PUBLIC    ::   rn_greenland_calving_fraction  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_antarctica_total_fw_flux   ! Constant total rate of freshwater input (kg/s) for Antarctica (if nn_coupled_iceshelf_fluxes=2)
+   REAL(wp), PUBLIC    ::   rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_iceshelf_fluxes_tolerance   ! Absolute tolerance for detecting differences in icesheet masses.
+	
    !!----------------------------------------------------------------------
    !!                     Surface atmospheric fields
    !!----------------------------------------------------------------------
@@ -193,6 +222,9 @@ CONTAINS
          &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) , cloud_fra(jpi,jpj),   &
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      &
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
+
+      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) ) 
+
       !
       ALLOCATE( e3t_m(jpi,jpj) , STAT=ierr(5) )
       !

src/OCE/SBC/sbc_phy.F90

@@ -864,7 +864,7 @@ CONTAINS
          zCh = zz0*ptst(ji,jj)/zdt
          zCe = zz0*pqst(ji,jj)/zdq
 
-         CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
+         CALL bulk_formula( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
             &              pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
             &              pTau(ji,jj), zQsen, zQlat )
 
@@ -879,7 +879,7 @@ CONTAINS
    END SUBROUTINE UPDATE_QNSOL_TAU
 
 
-   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
+   SUBROUTINE bulk_formula_sclr( pzu, pTs, pqs, pTa, pqa, &
       &                          pCd, pCh, pCe,           &
       &                          pwnd, pUb, ppa, prhoa,   &
       &                          pTau, pQsen, pQlat,      &
@@ -921,9 +921,9 @@ CONTAINS
 
       IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
 
-   END SUBROUTINE BULK_FORMULA_SCLR
+   END SUBROUTINE bulk_formula_sclr
 
-   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
+   SUBROUTINE bulk_formula_vctr( pzu, pTs, pqs, pTa, pqa, &
       &                          pCd, pCh, pCe,           &
       &                          pwnd, pUb, ppa, prhoa,   &
       &                          pTau, pQsen, pQlat,      &
@@ -957,7 +957,7 @@ CONTAINS
 
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
 
-         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
+         CALL bulk_formula_sclr( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
             &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  &
             &                    pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj),   &
             &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             &
@@ -966,7 +966,7 @@ CONTAINS
          IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
       END_2D
 
-   END SUBROUTINE BULK_FORMULA_VCTR
+   END SUBROUTINE bulk_formula_vctr
 
 
    FUNCTION alpha_sw_vctr( psst )

src/OCE/SBC/sbcblk.F90

@@ -117,8 +117,9 @@ MODULE sbcblk
    REAL(wp) ::   rn_stau_a      ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta
    REAL(wp) ::   rn_stau_b      !
    !
-   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation
+   REAL(dp)         ::   rn_pfac   ! multiplication factor for precipitation
    REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation
+   REAL(wp)         ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress                        !
    REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements
    REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements
    !
@@ -169,6 +170,7 @@ MODULE sbcblk
 
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "single_precision_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: sbcblk.F90 15551 2021-11-28 20:19:36Z gsamson $
@@ -219,7 +221,7 @@ CONTAINS
       INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable
       NAMELIST/namsbc_blk/ ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, ln_ANDREAS, &   ! bulk algorithm
          &                 rn_zqt, rn_zu, nn_iter_algo, ln_skin_cs, ln_skin_wl,       &
-         &                 rn_pfac, rn_efac,                                          &
+         &                 rn_pfac, rn_efac, rn_vfac,                                 &
          &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback
          &                 ln_humi_sph, ln_humi_dpt, ln_humi_rlh, ln_tair_pot,        &
          &                 ln_Cx_ice_cst, rn_Cd_i, rn_Ce_i, rn_Ch_i,                  &
@@ -413,6 +415,7 @@ CONTAINS
          WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
          WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
          WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
+         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
          WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
          WRITE(numout,*) '      use surface current feedback on wind stress         ln_crt_fbk   = ', ln_crt_fbk
          IF(ln_crt_fbk) THEN
@@ -575,7 +578,7 @@ CONTAINS
 
          CALL blk_oce_2(     theta_air_zt(:,:),                                    &   !   <<= in
             &                sf(jp_qlw  )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1),   &   !   <<= in
-            &                sf(jp_snow )%fnow(:,:,1), tsk_m,                      &   !   <<= in
+            &                CASTDP(sf(jp_snow )%fnow(:,:,1)), tsk_m,                      &   !   <<= in
             &                zsen, zlat, zevp )                                        !   <=> in out
       ENDIF
       !
@@ -654,9 +657,7 @@ CONTAINS
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) ::   zztmp                ! local variable
       REAL(wp) ::   zstmax, zstau
-#if defined key_cyclone
       REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
-#endif
       REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point
       REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
       REAL(wp), DIMENSION(jpi,jpj) ::   zcd_oce           ! momentum transfert coefficient over ocean
@@ -682,9 +683,10 @@ CONTAINS
       ! ----------------------------------------------------------------------------- !
 
       ! ... components ( U10m - U_oce ) at T-point (unmasked)
-#if defined key_cyclone
+
       zwnd_i(:,:) = 0._wp
       zwnd_j(:,:) = 0._wp
+#if defined key_cyclone
       CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj)
@@ -694,9 +696,15 @@ CONTAINS
       END_2D
 #else
       ! ... scalar wind module at T-point (not masked)
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
+      DO_2D( 0, 0, 0, 0 )
+         zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
+         zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
       END_2D
+      
+      CALL lbc_lnk( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
+      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
+      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
+         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
 #endif
       ! ----------------------------------------------------------------------------- !
       !      I   Solar FLUX                                                           !
@@ -808,7 +816,7 @@ CONTAINS
             rhoa(ji,jj) = rho_air( ztabs(ji,jj), q_zu(ji,jj), zpre(ji,jj) )
          END_2D
 
-         CALL BULK_FORMULA( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
+         CALL bulk_formula( rn_zu, zsspt(:,:), pssq(:,:), theta_zu(:,:), q_zu(:,:), &
             &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          &
             &               wndm(:,:), zU_zu(:,:), pslp(:,:), rhoa(:,:),       &
             &               taum(:,:), psen(:,:), plat(:,:),                   &
@@ -826,8 +834,13 @@ CONTAINS
                ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
                ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
 #else
-               ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
-               ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
+               IF ( rn_vfac > 0._wp ) THEN
+                   ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+                   ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+               ELSE
+                   ztau_i(ji,jj) = zztmp * pwndi(ji,jj)
+                   ztau_j(ji,jj) = zztmp * pwndj(ji,jj)
+               ENDIF
 #endif
             ELSE
                ztau_i(ji,jj) = 0._wp
@@ -868,11 +881,11 @@ CONTAINS
          CALL iom_put( "vtau_oce", ztau_j(:,:)*tmask(:,:,1) )  ! vtau at T-points!
 
          IF(sn_cfctl%l_prtctl) THEN
-            CALL prt_ctl( tab2d_1=pssq   , clinfo1=' blk_oce_1: pssq   : ', mask1=tmask )
-            CALL prt_ctl( tab2d_1=wndm   , clinfo1=' blk_oce_1: wndm   : ', mask1=tmask )
-            CALL prt_ctl( tab2d_1=utau   , clinfo1=' blk_oce_1: utau   : ', mask1=umask,   &
-               &          tab2d_2=vtau   , clinfo2='            vtau   : ', mask2=vmask )
-            CALL prt_ctl( tab2d_1=zcd_oce, clinfo1=' blk_oce_1: Cd     : ', mask1=tmask )
+            CALL prt_ctl( tab2d_1=CASTDP(pssq), clinfo1=' blk_oce_1: pssq   : ', mask1=tmask )
+            CALL prt_ctl( tab2d_1=CASTDP(wndm), clinfo1=' blk_oce_1: wndm   : ', mask1=tmask )
+            CALL prt_ctl( tab2d_1=CASTDP(utau), clinfo1=' blk_oce_1: utau   : ', mask1=umask,   &
+               &          tab2d_2=CASTDP(vtau), clinfo2='            vtau   : ', mask2=vmask )
+            CALL prt_ctl( tab2d_1=CASTDP(zcd_oce), clinfo1=' blk_oce_1: Cd     : ', mask1=tmask )
          ENDIF
          !
       ENDIF ! ln_blk / ln_abl
@@ -907,7 +920,7 @@ CONTAINS
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptair   ! potential temperature of air #LB: confirm!
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   pdqlw   ! downwelling longwave radiation at surface [W/m^2]
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   pprec
-      REAL(wp), INTENT(in), DIMENSION(:,:) ::   psnow
+      REAL(dp), INTENT(in), DIMENSION(:,:) ::   psnow
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptsk   ! SKIN surface temperature   [Celsius]
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   psen
       REAL(wp), INTENT(in), DIMENSION(:,:) ::   plat
@@ -969,11 +982,11 @@ CONTAINS
       ENDIF
       !
       IF(sn_cfctl%l_prtctl) THEN
-         CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw  : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=psen , clinfo1=' blk_oce_2: psen  : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=plat , clinfo1=' blk_oce_2: plat  : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=qns  , clinfo1=' blk_oce_2: qns   : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=emp  , clinfo1=' blk_oce_2: emp   : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(zqlw), clinfo1=' blk_oce_2: zqlw  : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(psen), clinfo1=' blk_oce_2: psen  : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(plat), clinfo1=' blk_oce_2: plat  : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(qns), clinfo1=' blk_oce_2: qns   : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(emp), clinfo1=' blk_oce_2: emp   : ', mask1=tmask )
       ENDIF
       !
    END SUBROUTINE blk_oce_2
@@ -1016,7 +1029,8 @@ CONTAINS
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
       REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature
-      REAL(wp) ::   zztmp1, zztmp2                ! temporary scalars
+      REAL(wp) ::   zwndi_t , zwndj_t 
+      REAL(wp) ::   zztmp0,zztmp1, zztmp2                ! temporary scalars
       REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zsipt ! temporary array
       !!---------------------------------------------------------------------
       !
@@ -1024,9 +1038,12 @@ CONTAINS
       !    Wind module relative to the moving ice ( U10m - U_ice )   !
       ! ------------------------------------------------------------ !
       ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
+      DO_2D( 0, 0, 0, 0 )
+         zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( puice(ji-1,jj  ) + puice(ji,jj) )  )
+         zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  )
+         wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t )* tmask(ji,jj,1)
       END_2D
+      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
       !
       ! potential sea-ice surface temperature [K]
       zsipt(:,:) = theta_exner( ptsui(:,:), pslp(:,:) )
@@ -1073,10 +1090,11 @@ CONTAINS
          !    Wind stress relative to nonmoving ice ( U10m )    !
          ! ---------------------------------------------------- !
          ! supress moving ice in wind stress computation as we don't know how to do it properly...
+         zztmp0 = rn_vfac * 0.5_wp
          DO_2D( 0, 1, 0, 1 )    ! at T point
             zztmp1        = rhoa(ji,jj) * Cd_ice(ji,jj) * wndm_ice(ji,jj)
-            putaui(ji,jj) = zztmp1 * pwndi(ji,jj)
-            pvtaui(ji,jj) = zztmp1 * pwndj(ji,jj)
+            putaui(ji,jj) = zztmp1 * ( pwndi(ji,jj) - zztmp0 * ( puice(ji-1,jj  ) + puice(ji,jj) ) )
+            pvtaui(ji,jj) = zztmp1 * ( pwndj(ji,jj) - zztmp0 * ( pvice(ji  ,jj-1) + pvice(ji,jj) ) )
          END_2D
 
          !#LB: saving the module, and x-y components, of the ai wind-stress at T-points: NOT weighted by the ice concentration !!!
@@ -1096,8 +1114,8 @@ CONTAINS
          END_2D
          CALL lbc_lnk( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
          !
-         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : ', mask1=umask   &
-            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ', mask2=vmask )
+         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=CASTDP(putaui), clinfo1=' blk_ice: putaui : ', mask1=umask   &
+            &                               , tab2d_2=CASTDP(pvtaui), clinfo2='          pvtaui : ', mask2=vmask )
       ELSE ! ln_abl
 
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
@@ -1109,7 +1127,7 @@ CONTAINS
 
       ENDIF ! ln_blk  / ln_abl
       !
-      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ', mask1=tmask )
+      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=CASTDP(wndm_ice), clinfo1=' blk_ice: wndm_ice : ', mask1=tmask )
       !
    END SUBROUTINE blk_ice_1
 
@@ -1310,18 +1328,18 @@ CONTAINS
          DO jl = 1, jpl
             zmsk(:,:,jl) = tmask(:,:,1)
          END DO
-         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', mask1=zmsk,   &
-            &         tab3d_2=z_qsb   , clinfo2=' z_qsb    : '         , mask2=zmsk, kdim=jpl)
-         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', mask1=zmsk,   &
-            &         tab3d_2=dqla_ice, clinfo2=' dqla_ice : '         , mask2=zmsk, kdim=jpl)
-         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', mask1=zmsk,   &
-            &         tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : '         , mask2=zmsk, kdim=jpl)
-         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', mask1=zmsk,   &
-            &         tab3d_2=qsr_ice , clinfo2=' qsr_ice  : '         , mask2=zmsk, kdim=jpl)
-         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', mask1=zmsk,   &
-            &         tab3d_2=qns_ice , clinfo2=' qns_ice  : '         , mask2=zmsk, kdim=jpl)
-         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', mask1=tmask,   &
-            &         tab2d_2=sprecip , clinfo2=' sprecip  : '         , mask2=tmask         )
+         CALL prt_ctl(tab3d_1=CASTDP(qla_ice) , clinfo1=' blk_ice: qla_ice  : ', mask1=zmsk,   &
+            &         tab3d_2=CASTDP(z_qsb)   , clinfo2=' z_qsb    : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=CASTDP(z_qlw)   , clinfo1=' blk_ice: z_qlw    : ', mask1=zmsk,   &
+            &         tab3d_2=CASTDP(dqla_ice), clinfo2=' dqla_ice : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=CASTDP(z_dqsb)  , clinfo1=' blk_ice: z_dqsb   : ', mask1=zmsk,   &
+            &         tab3d_2=CASTDP(z_dqlw)  , clinfo2=' z_dqlw   : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=CASTDP(dqns_ice), clinfo1=' blk_ice: dqns_ice : ', mask1=zmsk,   &
+            &         tab3d_2=CASTDP(qsr_ice) , clinfo2=' qsr_ice  : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab3d_1=CASTDP(ptsu)    , clinfo1=' blk_ice: ptsu     : ', mask1=zmsk,   &
+            &         tab3d_2=CASTDP(qns_ice) , clinfo2=' qns_ice  : '         , mask2=zmsk, kdim=jpl)
+         CALL prt_ctl(tab2d_1=CASTDP(tprecip) , clinfo1=' blk_ice: tprecip  : ', mask1=tmask,   &
+            &         tab2d_2=CASTDP(sprecip) , clinfo2=' sprecip  : '         , mask2=tmask         )
          DEALLOCATE(zmsk)
       ENDIF
 

src/OCE/SBC/sbcblk_algo_coare3p0.F90

@@ -169,9 +169,9 @@ CONTAINS
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
       !
       INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Cdn
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Chn
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Cen
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]

src/OCE/SBC/sbcblk_algo_coare3p6.F90

@@ -159,9 +159,9 @@ CONTAINS
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
       !
       INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Cdn
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Chn
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Cen
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]

src/OCE/SBC/sbcblk_algo_ecmwf.F90

@@ -165,9 +165,9 @@ CONTAINS
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
       !
       INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
-      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Cdn
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Chn
+      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Cen
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
       REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
@@ -417,27 +417,23 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
-      REAL(wp) :: zta, zx2, zx, ztmp, zpsi_unst, zpsi_stab, zstab, zc
+      REAL(wp) :: zta, zx2, zx, ztmp, zc
       !!----------------------------------------------------------------------------------
       zc = 5._wp/0.35_wp
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             !
             zta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!):
-
-            ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
-            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
-            zx  = SQRT(zx2)                          ! (1 - 16z)^0.25
-            ztmp = 1._wp + zx
-            zpsi_unst = LOG( 0.125_wp*ztmp*ztmp*(1._wp + zx2) ) - 2._wp*ATAN( zx ) + 0.5_wp*rpi
-
-            ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
-            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
-               &       - zta - 2._wp/3._wp*zc
             !
-            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
-            !
-            psi_m_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable
-               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
+            IF (zta < 0.0_wp) THEN
+               ! *** Unstable (Paulson 1970)    [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+               zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 - 16z)^0.5
+               zx  = SQRT(zx2)                        ! (1 - 16z)^0.25
+               ztmp = 1._wp + zx
+               psi_m_ecmwf(ji,jj) = LOG(0.125_wp*ztmp*ztmp*(1._wp + zx2)) - 2._wp*ATAN(zx) + 0.5_wp*rpi
+            ELSE
+               ! *** Stable                   [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+               psi_m_ecmwf(ji,jj) = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) - zta - 2._wp/3._wp*zc
+            ENDIF
             !
       END_2D
    END FUNCTION psi_m_ecmwf
@@ -458,7 +454,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
-      REAL(wp) ::  zta, zx2, zpsi_unst, zpsi_stab, zstab, zc
+      REAL(wp) ::  zta, zx2, zc
       !!----------------------------------------------------------------------------------
       zc = 5._wp/0.35_wp
       !
@@ -466,20 +462,15 @@ CONTAINS
             !
             zta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!):
             !
-            ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
-            zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
-            zpsi_unst = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
-            !
-            ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
-            zpsi_stab = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
+            IF (zta < 0.0_wp) THEN
+               ! *** Unstable (Paulson 1970)   [eq.3.20, Chap.3, p.33, IFS doc - Cy31r1] :
+               zx2 = SQRT( ABS(1._wp - 16._wp*zta) )  ! (1 -16z)^0.5
+               psi_h_ecmwf(ji,jj) = 2._wp*LOG( 0.5_wp*(1._wp + zx2) )
+            ELSE
+               ! *** Stable [eq.3.22, Chap.3, p.33, IFS doc - Cy31r1] :
+               psi_h_ecmwf(ji,jj) = -2._wp/3._wp*(zta - zc)*EXP(-0.35_wp*zta) &
                &       - ABS(1._wp + 2._wp/3._wp*zta)**1.5_wp - 2._wp/3._wp*zc + 1._wp
-            !
-            ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution...
-            !
-            zstab = 0.5_wp + SIGN(0.5_wp, zta) ! zta > 0 => zstab = 1
-            !
-            psi_h_ecmwf(ji,jj) =         zstab  * zpsi_stab &  ! (zta > 0) Stable
-               &              + (1._wp - zstab) * zpsi_unst    ! (zta < 0) Unstable
+            ENDIF
             !
       END_2D
    END FUNCTION psi_h_ecmwf

src/OCE/SBC/sbcclo.F90

@@ -159,7 +159,7 @@ MODULE sbcclo
          !
       END DO  ! jcs
 
-   END SUBROUTINE
+   END SUBROUTINE get_cssrcsurf
 
    SUBROUTINE get_cstrgsurf(kncs, kmaskcs, kmaskcsgrp, psurftrg, kcsgrp )
       !!-----------------------------------------------------------------------
@@ -211,7 +211,7 @@ MODULE sbcclo
          !
       END DO ! jcs
 
-   END SUBROUTINE
+   END SUBROUTINE get_cstrgsurf
 
    SUBROUTINE prt_csctl(kncs, psurfsrc, psurftrg, kcsgrp, cdcstype)
       !!-----------------------------------------------------------------------
@@ -245,7 +245,7 @@ MODULE sbcclo
          WRITE(numout,*)''
       END IF
 
-   END SUBROUTINE
+   END SUBROUTINE prt_csctl
 
    SUBROUTINE sbc_csupdate(kncs, kcsgrp, kmsk_src, kmsk_grp, psurfsrc, psurftrg, cdcstype, kmsk_opnsea, psurf_opnsea, pwcs, pqcs)
       !!-----------------------------------------------------------------------
@@ -319,7 +319,7 @@ MODULE sbcclo
          !
       END DO ! jcs
 
-   END SUBROUTINE
+   END SUBROUTINE sbc_csupdate
 
    SUBROUTINE alloc_csarr( klen, pvarsrc, pvartrg, kvargrp )
       !!-----------------------------------------------------------------------
@@ -347,6 +347,6 @@ MODULE sbcclo
       pvarsrc(:) = 0.e0_wp
       pvartrg(:) = 0.e0_wp
       kvargrp(:) = 0
-   END SUBROUTINE
+   END SUBROUTINE alloc_csarr
 
 END MODULE

src/OCE/SBC/sbccpl.F90

@@ -33,10 +33,18 @@ MODULE sbccpl
 #endif
    USE cpl_oasis3     ! OASIS3 coupling
    USE geo2ocean      !
-   USE oce     , ONLY : ts, uu, vv, ssh, fraqsr_1lev
+   USE oce   
+#if defined key_medusa
+   USE oce , ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl,  &
+                        PCO2a_in_cpl, Dust_in_cpl
+#endif
    USE ocealb         !
    USE eosbn2         !
    USE sbcrnf  , ONLY : l_rnfcpl
+   USE cpl_rnf_1d, ONLY: nn_cpl_river, cpl_rnf_1d_init, cpl_rnf_1d_to_2d   ! Variables used in 1D river outflow 
+#if defined key_medusa
+   USE par_trc , ONLY : ln_medusa
+#endif
 #if defined key_cice
    USE ice_domain_size, only: ncat
 #endif
@@ -48,6 +56,7 @@ MODULE sbccpl
    USE iom            ! NetCDF library
    USE lib_mpp        ! distribued memory computing library
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE timing
 
 #if defined key_oasis3
    USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut
@@ -129,8 +138,17 @@ MODULE sbccpl
    INTEGER, PARAMETER ::   jpr_icb    = 61
    INTEGER, PARAMETER ::   jpr_ts_ice = 62   ! Sea ice surface temp
    !!INTEGER, PARAMETER ::   jpr_qtrice = 63   ! Transmitted solar thru sea-ice
-
-   INTEGER, PARAMETER ::   jprcv      = 62   ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_grnm   = 63   ! Greenland ice mass 
+   INTEGER, PARAMETER ::   jpr_antm   = 64   ! Antarctic ice mass 
+   INTEGER, PARAMETER ::   jpr_rnf_1d = 65   ! 1D river runoff 
+   INTEGER, PARAMETER ::   jpr_qtr    = 66   ! Transmitted solar
+#if defined key_medusa
+   INTEGER, PARAMETER ::   jpr_atm_pco2 = 67 ! Incoming atm pCO2 flux
+   INTEGER, PARAMETER ::   jpr_atm_dust = 68 ! Incoming atm aggregate dust
+   INTEGER, PARAMETER ::   jprcv      = 69   ! total number of fields received
+#else
+   INTEGER, PARAMETER ::   jprcv      = 66   ! total number of fields received
+#endif
 
    INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
    INTEGER, PARAMETER ::   jps_toce   =  2   ! ocean temperature
@@ -170,17 +188,24 @@ MODULE sbccpl
    INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity
    INTEGER, PARAMETER ::   jps_sstfrz = 37   ! sea surface freezing temperature
    INTEGER, PARAMETER ::   jps_ttilyr = 38   ! sea ice top layer temp
-
-   INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent
-
-#if ! defined key_oasis3
-   ! Dummy variables to enable compilation when oasis3 is not being used
-   INTEGER                    ::   OASIS_Sent        = -1
-   INTEGER                    ::   OASIS_SentOut     = -1
-   INTEGER                    ::   OASIS_ToRest      = -1
-   INTEGER                    ::   OASIS_ToRestOut   = -1
+#if defined key_medusa
+   INTEGER, PARAMETER ::   jps_bio_co2 = 39    ! MEDUSA air-sea CO2 flux
+   INTEGER, PARAMETER ::   jps_bio_dms = 40    ! MEDUSA DMS surface concentration
+   INTEGER, PARAMETER ::   jps_bio_chloro = 41 ! MEDUSA chlorophyll surface concentration
+
+   INTEGER, PARAMETER ::   jpsnd      = 41   ! total number of fields sent 
+#else
+   INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent 
 #endif
 
+#if ! defined key_oasis3 
+   ! Dummy variables to enable compilation when oasis3 is not being used 
+   INTEGER                    ::   OASIS_Sent        = -1 
+   INTEGER                    ::   OASIS_SentOut     = -1 
+   INTEGER                    ::   OASIS_ToRest      = -1 
+   INTEGER                    ::   OASIS_ToRestOut   = -1 
+#endif 
+
    !                                  !!** namelist namsbc_cpl **
    TYPE ::   FLD_C                     !
       CHARACTER(len = 32) ::   cldes      ! desciption of the coupling strategy
@@ -193,20 +218,34 @@ MODULE sbccpl
    TYPE(FLD_C) ::   sn_snd_temp  , sn_snd_alb , sn_snd_thick, sn_snd_crt   , sn_snd_co2,  &
       &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
    !                                   ! Received from the atmosphere
-   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr,  &
+#if defined key_medusa
+   TYPE(FLD_C) ::   sn_snd_bio_co2, sn_snd_bio_dms, sn_snd_bio_chloro
+#endif
+   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr,  &
       &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf
-   !                                   ! Send to waves
-   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev
-   !                                   ! Received from waves
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf,      &
+                    sn_rcv_grnm, sn_rcv_antm
+#if defined key_medusa
+   TYPE(FLD_C) ::   sn_rcv_atm_pco2, sn_rcv_atm_dust
+#endif
+   ! Send to waves 
+   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
+
    TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, &
       &             sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd
+
+   ! Received from waves 
+   TYPE(FLD_C) ::   sn_rcv_tauwoc, sn_rcv_wfreq
+   ! Transmitted solar
+   TYPE(FLD_C) ::   sn_rcv_qtr
    !                                   ! Other namelist parameters
 !!   TYPE(FLD_C) ::   sn_rcv_qtrice
    INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
    LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models
                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
    LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
+   LOGICAL     ::   ln_couple_ocean_evap  ! Do we couple total (ocean+sea ice) evaporation (FALSE)
+                                          ! or ocean only evaporation (TRUE)
 
    TYPE ::   DYNARR
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3
@@ -215,14 +254,13 @@ MODULE sbccpl
    TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                ! all fields recieved from the atmosphere
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
-#if defined key_si3 || defined key_cice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time
-#endif
+
 
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) ::   nrcvinfo           ! OASIS info argument
 
    !! Substitution
 #  include "do_loop_substitute.h90"
+#  include "single_precision_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -274,22 +312,38 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
       !!
       NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             &
+         &                  ln_couple_ocean_evap,                                                      &
          &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  &
          &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  &
          &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  &
          &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  &
          &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_wstrf ,  &
          &                  sn_rcv_charn , sn_rcv_taw   , sn_rcv_bhd  , sn_rcv_tusd  , sn_rcv_tvsd,    &
+            sn_rcv_qtr   ,  &
          &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   &
          &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_icb  , sn_rcv_isf   , sn_rcv_ts_ice,  & !!, sn_rcv_qtrice
-         &                  sn_rcv_mslp
+         &                  sn_rcv_mslp, &
+                            sn_rcv_grnm  , sn_rcv_antm  ,                               &
+         &                  nn_coupled_iceshelf_fluxes  , ln_iceshelf_init_atmos ,                     &
+         &                  rn_greenland_total_fw_flux  , rn_greenland_calving_fraction  ,             &
+         &                  rn_antarctica_total_fw_flux , rn_antarctica_calving_fraction ,             &
+#if defined key_medusa
+         &                  rn_iceshelf_fluxes_tolerance,                                              &
+                                                                !  Add MEDUSA related fields to namelist
+                            sn_snd_bio_co2 , sn_snd_bio_dms, sn_snd_bio_chloro,                        &
+         &                  sn_rcv_atm_pco2, sn_rcv_atm_dust
+#else
+         &                  rn_iceshelf_fluxes_tolerance
+#endif
+
 
       !!---------------------------------------------------------------------
       !
       ! ================================ !
       !      Namelist informations       !
       ! ================================ !
-      !
+      ! 
+      IF (ln_timing) CALL timing_start('sbc_cpl_init')
       READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist' )
       !
@@ -306,6 +360,7 @@ CONTAINS
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
          WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux
+         WRITE(numout,*)'  ln_couple_ocean_evap                = ', ln_couple_ocean_evap
          WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
          WRITE(numout,*)'  received fields (mutiple ice categogies)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
@@ -320,10 +375,12 @@ CONTAINS
          WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
          WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
          WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
+         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')' 
+         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')' 
          WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')'
          WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
-!!       WRITE(numout,*)'      transmitted solar thru sea-ice  = ', TRIM(sn_rcv_qtrice%cldes), ' (', TRIM(sn_rcv_qtrice%clcat), ')'
+         WRITE(numout,*)'      transmitted solar               = ', TRIM(sn_rcv_qtr%cldes   ), ' (', TRIM(sn_rcv_qtr%clcat   ), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
          WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')'
          WRITE(numout,*)'      surface waves:'
@@ -346,6 +403,11 @@ CONTAINS
          WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
+#if defined key_medusa
+         WRITE(numout,*)'      bio co2 flux                    = ', TRIM(sn_snd_bio_co2%cldes),' (', TRIM(sn_snd_bio_co2%clcat), ')'
+         WRITE(numout,*)'      bio dms flux                    = ', TRIM(sn_snd_bio_dms%cldes),' (', TRIM(sn_snd_bio_dms%clcat), ')'
+         WRITE(numout,*)'      bio dms chlorophyll             = ', TRIM(sn_snd_bio_chloro%cldes), ' (', TRIM(sn_snd_bio_chloro%clcat), ')'
+#endif
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
          WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes  ), ' (', TRIM(sn_snd_cond%clcat  ), ')'
          WRITE(numout,*)'      meltponds fraction and depth    = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat  ), ')'
@@ -356,6 +418,17 @@ CONTAINS
          WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd
+#if defined key_medusa
+         WRITE(numout,*)'      atm pco2                        = ', TRIM(sn_rcv_atm_pco2%cldes),'(', TRIM(sn_rcv_atm_pco2%clcat   ), ')'
+         WRITE(numout,*)'      atm dust                        = ', TRIM(sn_rcv_atm_dust%cldes),'(', TRIM(sn_rcv_atm_dust%clcat),')'
+#endif
+          WRITE(numout,*)'  nn_coupled_iceshelf_fluxes          = ', nn_coupled_iceshelf_fluxes
+         WRITE(numout,*)'  ln_iceshelf_init_atmos              = ', ln_iceshelf_init_atmos
+         WRITE(numout,*)'  rn_greenland_total_fw_flux          = ', rn_greenland_total_fw_flux
+         WRITE(numout,*)'  rn_antarctica_total_fw_flux         = ', rn_antarctica_total_fw_flux
+         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction
+         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction
+         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance
       ENDIF
       IF( lwp .AND. ln_wave) THEN                        ! control print
       WRITE(numout,*)'      surface waves:'
@@ -390,7 +463,11 @@ CONTAINS
       !                 define the north fold type of lbc                  (srcv(:)%nsgn)
 
       ! default definitions of srcv
-      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
+      srcv(:)%laction = .FALSE. 
+      srcv(:)%clgrid = 'T' 
+      srcv(:)%nsgn = 1. 
+      srcv(:)%nct = 1 
+      srcv(:)%dimensions = 2 
 
       !                                                      ! ------------------------- !
       !                                                      ! ice and ocean wind stress !
@@ -444,7 +521,11 @@ CONTAINS
          srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
          srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
          srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
-         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
+         !srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
+ 	 ! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment
+         srcv(jpr_otx1)%laction = .TRUE.
+         srcv(jpr_oty1)%laction = .TRUE.
+         !
          srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
       CASE( 'T,I' )
          srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
@@ -509,15 +590,34 @@ CONTAINS
       !                                                      !     Runoffs & Calving     !
       !                                                      ! ------------------------- !
       srcv(jpr_rnf   )%clname = 'O_Runoff'
-      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
-         srcv(jpr_rnf)%laction = .TRUE.
+      srcv(jpr_rnf_1d   )%clname = 'ORunff1D' 
+      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' .OR. TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN  
+         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) srcv(jpr_rnf)%laction = .TRUE. 
+         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN 
+            srcv(jpr_rnf_1d)%laction = .TRUE. 
+            srcv(jpr_rnf_1d)%dimensions = 1 ! 1D field passed through coupler 
+         END IF 
          l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
          ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
       ENDIF
       !
-      srcv(jpr_cal)%clname = 'OCalving'   ;  IF( TRIM( sn_rcv_cal%cldes) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
+      srcv(jpr_cal   )%clname = 'OCalving'    
+      IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.      
+ 
+      srcv(jpr_grnm  )%clname = 'OGrnmass'  
+      IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' ) THEN 
+         srcv(jpr_grnm)%laction = .TRUE.  
+         srcv(jpr_grnm)%dimensions = 0 ! Scalar field
+      ENDIF 
+       
+      srcv(jpr_antm  )%clname = 'OAntmass' 
+      IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' ) THEN
+         srcv(jpr_antm)%laction = .TRUE.
+         srcv(jpr_antm)%dimensions = 0 ! Scalar field
+      ENDIF
+
       srcv(jpr_isf)%clname = 'OIcshelf'   ;  IF( TRIM( sn_rcv_isf%cldes) == 'coupled' )   srcv(jpr_isf)%laction = .TRUE.
       srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE.
 
@@ -593,6 +693,22 @@ CONTAINS
          IF(lwp) WRITE(numout,*) '   Atmospheric pco2 received from oasis '
          IF(lwp) WRITE(numout,*)
       ENDIF
+
+#if defined key_medusa   
+      !                                                      ! --------------------------------------- !
+      !                                                      ! Incoming CO2 and DUST fluxes for MEDUSA !
+      !                                                      ! --------------------------------------- !
+      srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
+      IF (TRIM(sn_rcv_atm_pco2%cldes) == 'medusa') THEN
+        srcv(jpr_atm_pco2)%laction = .TRUE.
+      END IF
+
+      srcv(jpr_atm_dust)%clname = 'OATMDUST'
+      IF (TRIM(sn_rcv_atm_dust%cldes) == 'medusa')  THEN
+        srcv(jpr_atm_dust)%laction = .TRUE.
+      END IF
+#endif
+
       !
       !                                                      ! ------------------------- !
       !                                                      ! Mean Sea Level Pressure   !
@@ -624,6 +740,21 @@ CONTAINS
 !!         ENDIF
 !!         srcv(jpr_qtrice)%laction = .TRUE.
 !!      ENDIF
+
+      !                                                      ! ------------------------- !
+      !                                                      !    transmitted solar      !   
+      !                                                      ! ------------------------- !
+      srcv(jpr_qtr )%clname = 'OQtr'
+      IF( TRIM(sn_rcv_qtr%cldes) == 'coupled' ) THEN
+         IF ( TRIM( sn_rcv_qtr%clcat ) == 'yes' ) THEN
+            srcv(jpr_qtr)%nct = nn_cats_cpl
+         ELSE
+            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtr%clcat should always be set to yes currently' )
+         ENDIF
+         srcv(jpr_qtr)%laction = .TRUE.
+      ENDIF
+
+
       !                                                      ! ------------------------- !
       !                                                      !    ice skin temperature   !
       !                                                      ! ------------------------- !
@@ -798,24 +929,6 @@ CONTAINS
          ENDIF
       ENDIF
 
-      ! =================================================== !
-      ! Allocate all parts of frcv used for received fields !
-      ! =================================================== !
-      DO jn = 1, jprcv
-         IF( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
-      END DO
-      ! Allocate taum part of frcv which is used even when not received as coupling field
-      IF( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
-      ! Allocate w10m part of frcv which is used even when not received as coupling field
-      IF( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
-      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
-      IF( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
-      IF( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
-      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
-      IF( k_ice /= 0 ) THEN
-         IF( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
-         IF( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
-      ENDIF
 
       ! ================================ !
       !     Define the send interface    !
@@ -825,8 +938,12 @@ CONTAINS
       !                 define the north fold type of lbc               (ssnd(:)%nsgn)
 
       ! default definitions of nsnd
-      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
-
+      ssnd(:)%laction = .FALSE. 
+      ssnd(:)%clgrid = 'T' 
+      ssnd(:)%nsgn = 1. 
+      ssnd(:)%nct = 1 
+      ssnd(:)%dimensions = 2 
+         
       !                                                      ! ------------------------- !
       !                                                      !    Surface temperature    !
       !                                                      ! ------------------------- !
@@ -977,6 +1094,25 @@ CONTAINS
       !                                                      ! ------------------------- !
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
       !
+
+#if defined key_medusa 
+      !                                                      ! ------------------------- !
+      !                                                      !   MEDUSA output fields    !
+      !                                                      ! ------------------------- !
+      ! Surface dimethyl sulphide from Medusa
+      ssnd(jps_bio_dms)%clname = 'OBioDMS'
+      IF( TRIM(sn_snd_bio_dms%cldes) == 'medusa' )    ssnd(jps_bio_dms )%laction = .TRUE.
+
+      ! Surface CO2 flux from Medusa
+      ssnd(jps_bio_co2)%clname = 'OBioCO2'
+      IF( TRIM(sn_snd_bio_co2%cldes) == 'medusa' )    ssnd(jps_bio_co2 )%laction = .TRUE.
+
+      ! Surface chlorophyll from Medusa
+      ssnd(jps_bio_chloro)%clname = 'OBioChlo'
+      IF( TRIM(sn_snd_bio_chloro%cldes) == 'medusa' ) ssnd(jps_bio_chloro )%laction = .TRUE.
+#endif
+
+
       !                                                      ! ------------------------- !
       !                                                      ! Sea surface freezing temp !
       !                                                      ! ------------------------- !
@@ -1104,10 +1240,52 @@ CONTAINS
          ENDIF
       ENDIF
 
+      ! Initialise 1D river outflow scheme 
+      nn_cpl_river = 1 
+      IF ( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) CALL cpl_rnf_1d_init   ! Coupled runoff using 1D array
+      
+      ! =================================================== !
+      ! Allocate all parts of frcv used for received fields !
+      ! =================================================== !
+      DO jn = 1, jprcv
+
+         IF ( srcv(jn)%laction ) THEN 
+            SELECT CASE( srcv(jn)%dimensions )
+            !
+            CASE( 0 )   ! Scalar field
+               ALLOCATE( frcv(jn)%z3(1,1,1) )
+               
+            CASE( 1 )   ! 1D field
+               ALLOCATE( frcv(jn)%z3(nn_cpl_river,1,1) )
+               
+            CASE DEFAULT  ! 2D (or pseudo 3D) field.
+               ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
+            END SELECT
+         END IF
+
+      END DO
+      ! Allocate taum part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
+      ! Allocate w10m part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
+      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
+      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
+      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
+      IF( k_ice /= 0 ) THEN
+         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
+         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
+      END IF
+
       !
       ! ================================ !
       !   initialisation of the coupler  !
       ! ================================ !
+
+      ! There's no point initialising the coupler if we've accumulated any errors in 
+      ! coupling field definitions or settings. 
+      IF (nstop > 0) CALL ctl_stop( 'STOP', 'sbc_cpl_init: Errors encountered in coupled field definitions' )
+
       CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
 
       IF(ln_usecplmask) THEN
@@ -1121,6 +1299,26 @@ CONTAINS
       ENDIF
       xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
       !
+      IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN 
+          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something 
+          ! more complicated could be done if required. 
+          greenland_icesheet_mask = 0.0 
+          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0 
+          antarctica_icesheet_mask = 0.0 
+          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0 
+  
+          IF( .not. ln_rstart ) THEN 
+             greenland_icesheet_mass = 0.0  
+             greenland_icesheet_mass_rate_of_change = 0.0  
+             greenland_icesheet_timelapsed = 0.0 
+             antarctica_icesheet_mass = 0.0  
+             antarctica_icesheet_mass_rate_of_change = 0.0  
+             antarctica_icesheet_timelapsed = 0.0 
+          ENDIF 
+ 
+      ENDIF 
+      !
+      IF (ln_timing) CALL timing_stop('sbc_cpl_init')
       !
    END SUBROUTINE sbc_cpl_init
 
@@ -1180,15 +1378,26 @@ CONTAINS
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
       INTEGER  ::   ji, jj, jn             ! dummy loop indices
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdt did not change since nit000)
+      INTEGER  ::   ikchoix
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztx2, zty2
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars
       REAL(wp) ::   zcoef                  ! temporary scalar
+      LOGICAL  ::   ll_wrtstp              ! write diagnostics?
       REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
+      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in 
+      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b 
+      REAL(wp) ::   zmask_sum, zepsilon    
       REAL(wp) ::   zzx, zzy               ! temporary variables
       REAL(wp) ::   r1_grau                ! = 1.e0 / (grav * rho0)
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
       !!----------------------------------------------------------------------
       !
+      !
+      IF (ln_timing) CALL timing_start('sbc_cpl_rcv')
+      !
+      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+            !
       IF( kt == nit000 ) THEN
       !   cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done
          ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
@@ -1208,7 +1417,16 @@ CONTAINS
       !                                                      ! ======================================================= !
       isec = ( kt - nit000 ) * NINT( rn_Dt )                      ! date of exchanges
       DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
-         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
+        IF( srcv(jn)%laction ) THEN  
+ 
+          IF ( srcv(jn)%dimensions <= 1 ) THEN 
+            CALL cpl_rcv_1d( jn, isec, frcv(jn)%z3, SIZE(frcv(jn)%z3), nrcvinfo(jn) ) 
+          ELSE 
+
+            CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) ) 
+          END IF 
+
+        END IF 
       END DO
 
       !                                                      ! ========================= !
@@ -1237,14 +1455,42 @@ CONTAINS
             !
             IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
                !                                                       ! (geographical to local grid -> rotate the components)
-               CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )
-               IF( srcv(jpr_otx2)%laction ) THEN
-                  CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )
+               IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN
+                  ! Temporary code for HadGEM3 - will be removed eventually.
+                  ! Only applies when we have only taux on U grid and tauy on V grid
+
+                  !RSRH these MUST be initialised because the halos are not explicitly set 
+                  ! but they're passed to repcmo and used directly in calculations, so if 
+                  ! they point at junk in memory then bad things will happen!
+                  ! (You can prove this by running with preset NaNs). 
+                  ztx(:,:)=0.0
+                  zty(:,:)=0.0
+
+                  DO_2D( 0, 0, 0, 0 )                                       
+                          ztx(ji,jj)=0.25*vmask(ji,jj,1)                &
+                             *(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1)    &
+                             +frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))
+                          zty(ji,jj)=0.25*umask(ji,jj,1)                &
+                             *(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1)    &
+                             +frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))
+                  END_2D
+	                       
+                  ikchoix = 1
+                  CALL repcmo (frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)
+                  CALL lbc_lnk ('jpr_otx1', ztx2,'U', -1. )
+                  CALL lbc_lnk ('jpr_oty1', zty2,'V', -1. )
+                  frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)
+                  frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)
                ELSE
-                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
+                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
+                  frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
+                  IF( srcv(jpr_otx2)%laction ) THEN
+                     CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
+                  ELSE
+                     CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty )
+                  ENDIF
+                  frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid 
                ENDIF
-               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
-               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid
             ENDIF
             !
             IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
@@ -1334,6 +1580,12 @@ CONTAINS
          !
       ENDIF
 
+#if defined key_medusa
+      IF (ln_medusa) THEN
+        IF( srcv(jpr_atm_pco2)%laction) PCO2a_in_cpl(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
+        IF( srcv(jpr_atm_dust)%laction) Dust_in_cpl(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
+      ENDIF
+#endif
       !                                                      ! ================== !
       !                                                      ! atmosph. CO2 (ppm) !
       !                                                      ! ================== !
@@ -1538,8 +1790,10 @@ CONTAINS
          !
          IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus ! remove heat content associated to iceberg melting
          !
-         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
-         ELSE                   ;   qns(:,:) =                              zqns(:,:)
+         IF( ln_mixcpl ) THEN   
+            qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
+         ELSE                   
+            qns(:,:) =                              zqns(:,:)
          ENDIF
 
          !                                                       ! solar flux over the ocean          (qsr)
@@ -1558,7 +1812,78 @@ CONTAINS
          IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
          !
       ENDIF
+
+      zepsilon = rn_iceshelf_fluxes_tolerance
+
+      IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
+     
+         ! This is a zero dimensional, single value field.
+         zgreenland_icesheet_mass_in =  frcv(jpr_grnm)%z3(1,1,1)
+           
+         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
+         IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
+            ! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
+            ! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
+            zgreenland_icesheet_mass_b = zgreenland_icesheet_mass_in
+            greenland_icesheet_mass = zgreenland_icesheet_mass_in
+         ENDIF
+         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
+            zgreenland_icesheet_mass_b = greenland_icesheet_mass
+           
+            ! Only update the mass if it has increased.
+            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
+               greenland_icesheet_mass = zgreenland_icesheet_mass_in
+            ENDIF
+           
+            IF( zgreenland_icesheet_mass_b /= 0.0 ) &
+           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed 
+            greenland_icesheet_timelapsed = 0.0_wp       
+         ENDIF
+         IF(lwp .AND. ll_wrtstp) THEN
+            WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
+            WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass
+            WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
+            WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
+         ENDIF
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         greenland_icesheet_mass_rate_of_change = rn_greenland_total_fw_flux
+      ENDIF
+      !                                                        ! land ice masses : Antarctica
+      IF( srcv(jpr_antm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
+         
+         ! This is a zero dimensional, single value field.
+         zantarctica_icesheet_mass_in = frcv(jpr_antm)%z3(1,1,1)
+           
+         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
+         IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
+            ! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
+            ! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
+            zantarctica_icesheet_mass_b = zantarctica_icesheet_mass_in
+            antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+         ENDIF
+         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
+            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
+           
+            ! Only update the mass if it has increased.
+            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
+               antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+            END IF
+           
+            IF( zantarctica_icesheet_mass_b /= 0.0 ) &
+          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed 
+            antarctica_icesheet_timelapsed = 0.0_wp       
+         ENDIF
+         IF(lwp .AND. ll_wrtstp) THEN
+            WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
+            WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass
+            WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
+            WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
+         ENDIF
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
+      ENDIF
       !
+      IF (ln_timing) CALL timing_stop('sbc_cpl_rcv')
    END SUBROUTINE sbc_cpl_rcv
 
 
@@ -1672,7 +1997,7 @@ CONTAINS
                p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
                p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
             END_2D
-            CALL lbc_lnk( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
+            CALL lbc_lnk( 'sbccpl', p_taui, 'U',  -1._wp, p_tauj, 'V',  -1._wp )
          END SELECT
 
       ENDIF
@@ -1758,17 +2083,32 @@ CONTAINS
       !!----------------------------------------------------------------------
       !
 #if defined key_si3 || defined key_cice
-      !
-      IF( kt == nit000 ) THEN
-         ! allocate ice fractions from last coupling time here and not in sbc_cpl_init because of jpl
-         IF( .NOT.ALLOCATED(a_i_last_couple) )   ALLOCATE( a_i_last_couple(jpi,jpj,jpl) )
-         ! initialize to a_i for the 1st time step
-         a_i_last_couple(:,:,:) = a_i(:,:,:)
-      ENDIF
+
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
       ziceld(:,:) = 1._wp - picefr(:,:)
       zcptn (:,:) = rcp * sst_m(:,:)
+
+
+      !                                                      ! ========================= !
+      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
+      !                                                      ! ========================= !
+      CASE ('coupled')
+         IF (ln_scale_ice_flux) THEN
+            WHERE( a_i(:,:,:) > 1.e-10_wp )
+               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+            ELSEWHERE
+               qml_ice(:,:,:) = 0.0_wp
+               qcn_ice(:,:,:) = 0.0_wp
+            END WHERE
+         ELSE
+            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
+            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
+         ENDIF
+      END SELECT
+
+
       !
       !                                                      ! ========================= !
       !                                                      !    freshwater budget      !   (emp_tot)
@@ -1782,7 +2122,9 @@ CONTAINS
       CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
          zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
          ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
-         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
+         IF (.not. ln_couple_ocean_evap ) THEN
+            zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
+         END IF
       CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
          zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
@@ -1790,7 +2132,9 @@ CONTAINS
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
       CASE( 'none'      )       ! Not available as for now: needs additional coding below when computing zevap_oce
       !                         ! since fields received are not defined with none option
-         CALL ctl_stop('STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl')
+         CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
+      CASE default              ! Default
+         CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
       END SELECT
 
       ! --- evaporation over ice (kg/m2/s) --- !
@@ -1836,19 +2180,30 @@ CONTAINS
 
       ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
       zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( picefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
-      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
 
-      ! --- evaporation over ocean (used later for qemp) --- !
-      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
+      IF ( ln_couple_ocean_evap ) THEN
+         zemp_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_rain)%z3(:,:,1) & !Ocean evap minus rain (as all rain goes straight to ocean in GC5)
+                        & - zsprecip(:,:) * ( 1._wp - zsnw(:,:) )        !subtract snow in leads after correction for blowing snow
+         zemp_tot(:,:) = zemp_oce(:,:) + zemp_ice(:,:)
+         zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1)
+      ELSE
+         zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
+
+         ! --- evaporation over ocean (used later for qemp) --- !
+         zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
+      END IF
 
       ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
       ! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
       zdevap_ice(:,:) = 0._wp
 
       ! --- Continental fluxes --- !
-      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
+      IF( srcv(jpr_rnf)%laction ) THEN   ! 2D runoffs (included in emp later on)
          rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
       ENDIF
+      IF( srcv(jpr_rnf_1d)%laction ) THEN ! 1D runoff
+         CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:))
+      ENDIF
       IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot and emp_oce)
          zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
          zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
@@ -1928,23 +2283,7 @@ CONTAINS
 !!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
 !!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
       !
-      !                                                      ! ========================= !
-      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
-      !                                                      ! ========================= !
-      CASE ('coupled')
-         IF (ln_scale_ice_flux) THEN
-            WHERE( a_i(:,:,:) > 1.e-10_wp )
-               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
-               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
-            ELSEWHERE
-               qml_ice(:,:,:) = 0.0_wp
-               qcn_ice(:,:,:) = 0.0_wp
-            END WHERE
-         ELSE
-            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
-            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
-         ENDIF
-      END SELECT
+
       !
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
@@ -2008,8 +2347,10 @@ CONTAINS
       ! --- non solar flux over ocean --- !
       !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
       zqns_oce = 0._wp
-      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
 
+
+
+      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
       ! Heat content per unit mass of snow (J/kg)
       WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
       ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:)
@@ -2022,7 +2363,12 @@ CONTAINS
 
       ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
       DO jl = 1, jpl
-         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
+         ! zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
+
+         ! How much enthalpy is stored in sublimating snow and ice. At this stage we don't know if it is snow or ice that is
+         ! sublimating so we will use the combined snow and ice layer temperature t1_ice.
+         zqevap_ice(:,:,jl) = -zevap_ice(:,:,jl) * ( ( rt0 - t1_ice(:,:,jl) ) * rcpi + rLfus )
+
       END DO
 
       ! --- heat flux associated with emp (W/m2) --- !
@@ -2040,6 +2386,7 @@ CONTAINS
       IF( ln_mixcpl ) THEN
          qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
          qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
+
          DO jl=1,jpl
             qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
             qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
@@ -2084,9 +2431,14 @@ CONTAINS
       IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * rLfus ) ! latent heat from icebergs melting
       IF (        iom_use('hflx_rain_cea') )    &                                                    ! heat flux from rain (cell average)
          &   CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
-      IF (        iom_use('hflx_evap_cea') )    &                                                    ! heat flux from evap (cell average)
-         &   CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) )  &
-         &                                  * zcptn(:,:) * tmask(:,:,1) )
+
+      IF ( ln_couple_ocean_evap ) THEN
+         IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , frcv(jpr_tevp)%z3(:,:,1) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
+      ELSE
+         IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) )  &
+           &                         * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
+      END IF
+
       IF (        iom_use('hflx_prec_cea') )    &                                                    ! heat flux from all precip (cell avg)
          &   CALL iom_put('hflx_prec_cea' ,    sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  &
          &                                 + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
@@ -2213,29 +2565,30 @@ CONTAINS
          !
       ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
          !
-!!         SELECT CASE( TRIM( sn_rcv_qtrice%cldes ) )
-!!            !
-!!            !      ! ===> here we receive the qtr_ice_top array from the coupler
-!!         CASE ('coupled')
-!!            IF (ln_scale_ice_flux) THEN
-!!               WHERE( a_i(:,:,:) > 1.e-10_wp )
-!!                  zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
-!!               ELSEWHERE
-!!                  zqtr_ice_top(:,:,:) = 0.0_wp
-!!               ENDWHERE
-!!            ELSE
-!!               zqtr_ice_top(:,:,:) = frcv(jpr_qtrice)%z3(:,:,:)
-!!            ENDIF
-!!           
-!!            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
-!!            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
-!!            zqsr_tot(:,:)   = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
-!!            
-!!            !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
-!!         CASE ('none')
-         zqtr_ice_top(:,:,:) = 0._wp
-!!         END SELECT
-            !
+         SELECT CASE( TRIM( sn_rcv_qtr%cldes ) )
+         !
+         !      ! ===> here we receive the qtr_ice_top array from the coupler
+         CASE ('coupled')
+
+            IF (ln_scale_ice_flux) THEN
+               WHERE( a_i(:,:,:) > 0.0_wp ) zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:)) * a_i_last_couple(:,:,:) / a_i(:,:,:)
+               WHERE( a_i(:,:,:) <= 0.0_wp ) zqtr_ice_top(:,:,:) = 0.0_wp
+            ELSE
+               zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:))
+            ENDIF
+
+            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
+            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
+            zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
+
+         !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
+         CASE ('none')
+            zqtr_ice_top(:,:,:) = 0._wp
+         CASE default
+            CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_qtr%cldes' )
+       END SELECT
+
+         !
       ENDIF
 
       IF( ln_mixcpl ) THEN
@@ -2259,6 +2612,44 @@ CONTAINS
       IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
       ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
 
+      IF( ln_mixcpl ) THEN
+         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
+         DO jl = 1, jpl
+            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
+         END DO
+      ELSE
+         qsr_tot(:,:  ) = zqsr_tot(:,:  )
+         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
+      ENDIF
+
+      !                                                      ! ========================= !
+      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
+      !                                                      ! ========================= !
+      CASE ('coupled')
+         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
+            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
+         ELSE
+            ! Set all category values equal for the moment
+            DO jl=1,jpl
+               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
+            ENDDO
+         ENDIF
+      CASE( 'none' ) 
+         zdqns_ice(:,:,:) = 0._wp
+      CASE default
+         CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_dqnsdt%cldes' )
+      END SELECT
+      
+      IF( ln_mixcpl ) THEN
+         DO jl=1,jpl
+            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
+         ENDDO
+      ELSE
+         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
+      ENDIF
+
+
       !                                                      ! ================== !
       !                                                      !   ice skin temp.   !
       !                                                      ! ================== !
@@ -2296,13 +2687,18 @@ CONTAINS
       INTEGER, INTENT(in) ::   kt
       INTEGER, INTENT(in) ::   Kbb, Kmm    ! ocean model time level index
       !
-      INTEGER ::   ji, jj, jl   ! dummy loop indices
+      INTEGER ::   ji, jj, jl   ! dummy loop indices      
+      INTEGER ::   ikchoix
       INTEGER ::   isec, info   ! local integer
       REAL(wp) ::   zumax, zvmax
       REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp5, ztmp6 ! RSRH temporary work arrays 
+                                                         ! to avoid intent conflicts in repcmo calls
       !!----------------------------------------------------------------------
       !
+      IF (ln_timing) CALL timing_start('sbc_cpl_snd')
+
       isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges
       info = OASIS_idle
 
@@ -2316,12 +2712,15 @@ CONTAINS
             ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
          ELSE
             ! we must send the surface potential temperature
-            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( ts(:,:,1,jp_tem,Kmm), ts(:,:,1,jp_sal,Kmm) )
-            ELSE                   ;   ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
+            IF( l_useCT )  THEN    
+                ztmp1(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
+            ELSE                   
+                ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
             ENDIF
             !
             SELECT CASE( sn_snd_temp%cldes)
-            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
+            CASE( 'oce only'             )   
+                ztmp1(:,:) =   ztmp1(:,:) + rt0
             CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
                SELECT CASE( sn_snd_temp%clcat )
                CASE( 'yes' )
@@ -2563,18 +2962,37 @@ CONTAINS
       !                                                      ! ------------------------- !
       !                                                      !  CO2 flux from PISCES     !
       !                                                      ! ------------------------- !
-      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN
-         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s
-         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info )
+      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN  
+         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s 
+         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info ) 
+      ENDIF 
+#if defined key_medusa 
+      !
+      IF (ln_medusa) THEN
+      !                                                      ! ---------------------------------------------- !
+      !                                                      !  CO2 flux, DMS and chlorophyll from MEDUSA     !
+      !                                                      ! ---------------------------------------------- !
+         IF ( ssnd(jps_bio_co2)%laction ) THEN
+            CALL cpl_snd( jps_bio_co2, isec, RESHAPE( CO2Flux_out_cpl, (/jpi,jpj,1/) ), info )
+         ENDIF
+
+         IF ( ssnd(jps_bio_dms)%laction )  THEN
+            CALL cpl_snd( jps_bio_dms, isec, RESHAPE( DMS_out_cpl, (/jpi,jpj,1/) ), info )
+         ENDIF
+
+         IF ( ssnd(jps_bio_chloro)%laction )  THEN
+            CALL cpl_snd( jps_bio_chloro, isec, RESHAPE( chloro_out_cpl, (/jpi,jpj,1/) ), info )
+         ENDIF
       ENDIF
       !
+#endif
       !                                                      ! ------------------------- !
       IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
          !                                                   ! ------------------------- !
          !
          !                                                  j+1   j     -----V---F
          ! surface velocity always sent from T point                     !       |
-         !                                                        j      |   T   U
+         !  [except for HadGEM3]                                  j      |   T   U
          !                                                               |       |
          !                                                   j    j-1   -I-------|
          !                                               (for I)         |       |
@@ -2586,11 +3004,19 @@ CONTAINS
          ELSE
             SELECT CASE( TRIM( sn_snd_crt%cldes ) )
             CASE( 'oce only'             )      ! C-grid ==> T
-               DO_2D( 0, 0, 0, 0 )
-                  zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )
-                  zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji  ,jj-1,1,Kmm) )
-               END_2D
-            CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T
+               IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
+                  DO_2D( 0, 0, 0, 0 )
+                        zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )
+                        zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji  ,jj-1,1,Kmm) )
+                  END_2D
+               ELSE
+! Temporarily Changed for UKV
+                  DO_2D( 0, 0, 0, 0 )
+                        zotx1(ji,jj) = uu(ji,jj,1,Kmm)
+                        zoty1(ji,jj) = vv(ji,jj,1,Kmm)
+                  END_2D
+               ENDIF
+           CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T
                DO_2D( 0, 0, 0, 0 )
                   zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)
                   zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)
@@ -2606,6 +3032,7 @@ CONTAINS
                      &         + 0.5 * ( v_ice(ji,jj  )     + v_ice(ji  ,jj-1  )     ) *  fr_i(ji,jj)
                END_2D
             END SELECT
+               
             CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
             !
          ENDIF
@@ -2613,15 +3040,44 @@ CONTAINS
          !
          IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
             !                                                                     ! Ocean component
-            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
-            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
-            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
-            zoty1(:,:) = ztmp2(:,:)
-            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
-               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
-               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
-               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
-               zity1(:,:) = ztmp2(:,:)
+            IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
+               CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
+               CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
+               zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
+               zoty1(:,:) = ztmp2(:,:)
+               IF( ssnd(jps_ivx1)%laction ) THEN                                  ! Ice component
+                  CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
+                  CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
+                  zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
+                  zity1(:,:) = ztmp2(:,:)
+               ENDIF
+            ELSE
+               ! Temporary code for HadGEM3 - will be removed eventually.
+               ! Only applies when we want uvel on U grid and vvel on V grid
+               ! Rotate U and V onto geographic grid before sending.
+	
+              DO_2D( 0, 0, 0, 0 )
+                     ztmp1(ji,jj)=0.25*vmask(ji,jj,1)                  &
+                          *(zotx1(ji,jj)+zotx1(ji-1,jj)    &
+                          +zotx1(ji,jj+1)+zotx1(ji-1,jj+1))
+                     ztmp2(ji,jj)=0.25*umask(ji,jj,1)                  &
+                          *(zoty1(ji,jj)+zoty1(ji+1,jj)    &
+                          +zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
+               END_2D
+              	               
+               ikchoix = -1
+               ! zotx1 and zoty1 are input only to repcmo while ztmp5 and ztmp6
+               ! are the newly calculated (output) values.
+               ! Don't make the mistake of using zotx1 and zoty1 twice in this
+               ! call for both input and output fields since it creates INTENT
+               ! conflicts. 
+               CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,ztmp5,ztmp6,ikchoix)
+               zotx1(:,:)=ztmp5(:,:)
+               zoty1(:,:)=ztmp6(:,:)
+
+               ! Ensure any N fold and wrap columns are updated. 
+               CALL lbc_lnk( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp )
+
             ENDIF
          ENDIF
          !
@@ -2750,11 +3206,11 @@ CONTAINS
       ENDIF
       !                                                        ! SSS
       IF( ssnd(jps_soce  )%laction )  THEN
-         CALL cpl_snd( jps_soce  , isec, RESHAPE ( ts(:,:,1,jp_sal,Kmm), (/jpi,jpj,1/) ), info )
+         CALL cpl_snd( jps_soce  , isec, CASTSP(RESHAPE ( ts(:,:,1,jp_sal,Kmm), (/jpi,jpj,1/) )), info )
       ENDIF
       !                                                        ! first T level thickness
       IF( ssnd(jps_e3t1st )%laction )  THEN
-         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t(:,:,1,Kmm)   , (/jpi,jpj,1/) ), info )
+ CALL cpl_snd( jps_e3t1st, isec, CASTSP(RESHAPE ( e3t(:,:,1,Kmm) , (/jpi,jpj,1/) )), info )
       ENDIF
       !                                                        ! Qsr fraction
       IF( ssnd(jps_fraqsr)%laction )  THEN
@@ -2782,6 +3238,7 @@ CONTAINS
       IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info)
 #endif
       !
+      IF (ln_timing) CALL timing_stop('sbc_cpl_snd')
    END SUBROUTINE sbc_cpl_snd
 
    !!======================================================================

src/OCE/SBC/sbcfwb.F90

@@ -21,6 +21,9 @@ MODULE sbcfwb
    USE phycst         ! physical constants
    USE sbcrnf         ! ocean runoffs
    USE sbcssr         ! Sea-Surface damping terms
+#if defined key_agrif
+   USE agrif_oce , ONLY: Kmm_a
+#endif
    !
    USE in_out_manager ! I/O manager
    USE iom            ! IOM
@@ -37,8 +40,13 @@ MODULE sbcfwb
    REAL(wp) ::   rn_fwb0   ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only)
    REAL(wp) ::   a_fwb     ! annual domain averaged freshwater budget from the previous year
    REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget from the year before or at initial state
-   REAL(wp) ::   a_fwb_ini ! initial domain averaged freshwater budget
    REAL(wp) ::   area      ! global mean ocean surface (interior domain)
+   REAL(wp) ::   emp_corr  ! current, globally constant, emp correction
+#if defined key_agrif
+!$AGRIF_DO_NOT_TREAT
+   REAL(wp) ::   agrif_sum ! global sum used for agrif
+!$AGRIF_END_DO_NOT_TREAT
+#endif
 
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -67,8 +75,9 @@ CONTAINS
       INTEGER, INTENT( in ) ::   Kmm      ! ocean time level index
       !
       INTEGER  ::   ios, inum, ikty       ! local integers
-      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp                ! local scalars
-      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread, zcoef          !   -      -
+      INTEGER  ::   ji, jj, istart, iend, jstart, jend
+      REAL(wp) ::   z_fwf, z_fwf_nsrf, zsum_fwf, zsum_erp         ! local scalars
+      REAL(wp) ::   zsurf_neg, zsurf_pos, zsurf_tospread          !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_neg, ztmsk_pos, z_wgt ! 2D workspaces
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      -
       REAL(wp)   ,DIMENSION(1) ::   z_fwfprv  
@@ -99,8 +108,20 @@ CONTAINS
          !
          IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
          IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' )
+#if defined key_agrif
+         IF((kn_fwb == 3).AND.(Agrif_maxlevel()/=0)) CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 not yet implemented with AGRIF zooms ' ) 
+#endif
          !
-         area = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface
+         IF ( Agrif_Root() ) THEN
+#if defined key_agrif
+            agrif_sum = 0.0_wp                         ! set work scalar to 0
+            CALL Agrif_step_child(glob_sum_area_agrif) ! integrate over grids 
+            area = agrif_sum
+#else
+            area = glob_sum( 'sbcfwb', e1e2t(:,:) * tmask(:,:,1))         ! interior global domain surface
+#endif
+            IF (lwp) WRITE(numout,*) 'Domain area (km**2):', area/1000._wp/1000._wp
+         ENDIF
          ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes
          ! and in case of no melt, it can generate HSSW.
          !
@@ -116,77 +137,173 @@ CONTAINS
       !
       CASE ( 1 )                             !==  global mean fwf set to zero  ==!
          !
+#if ! defined key_agrif
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            y_fwfnow(1) = local_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
-            CALL mpp_delay_sum( 'sbcfwb', 'fwb', y_fwfnow(:), z_fwfprv(:), kt == nitend - nn_fsbc + 1 )
-            z_fwfprv(1) = z_fwfprv(1) / area
-            zcoef = z_fwfprv(1) * rcp
-            emp(:,:) = emp(:,:) - z_fwfprv(1)        * tmask(:,:,1)
-            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
-            ! outputs
-            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', zcoef * sst_m(:,:) * tmask(:,:,1) )
-            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1)        * tmask(:,:,1) )
+            z_fwfprv(1) = glob_sum( 'sbcfwb', e1e2t(:,:) * (  emp(:,:) - rnf(:,:)               & 
+                        &                                   - fwfisf_cav(:,:) - fwfisf_par(:,:) & 
+                        &                                   - snwice_fmass(:,:) ) ) / area
+            emp_corr = -z_fwfprv(1)
+            emp(:,:) = emp(:,:) + emp_corr * tmask(:,:,1)
+            qns(:,:) = qns(:,:) - emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+         ENDIF
+#else                                            
+         IF ( Agrif_Root() ) THEN
+            IF ( Agrif_maxlevel()==0 ) THEN
+               ! No child grid, correct "now" fluxes (i.e. as in the "no agrif" case)
+               IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+                  emp_corr = -glob_sum( 'sbcfwb', e1e2t(:,:) * (  emp(:,:) - rnf(:,:)                & 
+                             &                                   - fwfisf_cav(:,:) - fwfisf_par(:,:) & 
+                             &                                   - snwice_fmass(:,:) ) ) / area 
+              ENDIF
+            ELSE
+               !
+               ! Volume is here corrected according to the budget computed in the past, e.g. between
+               ! the last two consecutive calls to the surface module. Hence, the volume is allowed to drift slightly during
+               ! the current time step. 
+               !
+               IF( kt == nit000 ) THEN                                                               ! initialisation
+                  !                                                                                  ! 
+                  IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb',      ldstop = .FALSE. ) > 0     & ! read from restart file
+                     &           .AND. iom_varid( numror, 'emp_corr',   ldstop = .FALSE. ) > 0 ) THEN
+                     IF(lwp)   WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file'
+                     CALL iom_get( numror, 'a_fwb'     , a_fwb    )
+                     CALL iom_get( numror, 'emp_corr'  , emp_corr )
+                     !                                                                                 
+                  ELSE                                                                               
+                     emp_corr = 0._wp
+                     a_fwb    = 999._wp
+                     a_fwb_b  = 999._wp   
+                  END IF
+                  !
+                  IF(lwp)   WRITE(numout,*)
+                  IF(lwp)   WRITE(numout,*)'sbc_fwb : initial freshwater correction flux = ', emp_corr , 'kg/m2/s'
+                  !
+               ENDIF
+               !
+               IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+                  a_fwb_b   = a_fwb                            ! time swap
+                  agrif_sum = 0._wp                            ! set work scalar to 0
+                  CALL Agrif_step_child(glob_sum_volume_agrif) ! integrate over grids 
+                  a_fwb = agrif_sum * rho0 / area
+                  IF ( a_fwb_b == 999._wp ) a_fwb_b = a_fwb
+                  emp_corr = (a_fwb - a_fwb_b) / ( rn_Dt * REAL(kn_fsbc, wp) ) + emp_corr
+               ENDIF    
+               !   
+            ENDIF 
+         ELSE ! child grid if any
+            IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+                emp_corr  =  Agrif_parent(emp_corr) 
+            ENDIF  
          ENDIF
          !
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN        ! correct the freshwater fluxes on all grids
+            emp(:,:) = emp(:,:) + emp_corr * tmask(:,:,1)
+            qns(:,:) = qns(:,:) - emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+         ENDIF
+
+         IF ( Agrif_Root() ) THEN
+            ! Output restart information (root grid only)
+            IF( lrst_oce ) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
+               IF(lwp) WRITE(numout,*) '~~~~'
+               CALL iom_rstput( kt, nitrst, numrow, 'a_fwb'  ,   a_fwb   )
+               CALL iom_rstput( kt, nitrst, numrow, 'emp_corr', emp_corr )
+            END IF
+            !
+            IF( kt == nitend .AND. lwp ) THEN
+               WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', emp_corr  , 'kg/m2/s'
+            END IF
+         END IF
+#endif
+         ! outputs
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -emp_corr * tmask(:,:,1) )
+         ENDIF   
+         !
       CASE ( 2 )                             !==  fw adjustment based on fw budget at the end of the previous year  ==!
          !                                                simulation is supposed to start 1st of January
-         IF( kt == nit000 ) THEN                                                                 ! initialisation
-            !                                                                                    ! set the fw adjustment (a_fwb)
-            IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0     &      !    as read from restart file
-               &           .AND. iom_varid( numror, 'a_fwb',   ldstop = .FALSE. ) > 0 ) THEN
-               IF(lwp)   WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file'
-               CALL iom_get( numror, 'a_fwb_b', a_fwb_b )
-               CALL iom_get( numror, 'a_fwb'  , a_fwb )
+         IF ( Agrif_Root() ) THEN
+            IF( kt == nit000 ) THEN                                                                 ! initialisation
+               !                                                                                    ! 
+               IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb',      ldstop = .FALSE. ) > 0     &   ! read from restart file
+                  &           .AND. iom_varid( numror, 'a_fwb_b',    ldstop = .FALSE. ) > 0     & 
+                  &           .AND. iom_varid( numror, 'emp_corr',   ldstop = .FALSE. ) > 0 ) THEN
+                  IF(lwp)   WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file'
+                  CALL iom_get( numror, 'a_fwb'     , a_fwb    )
+                  CALL iom_get( numror, 'a_fwb_b'   , a_fwb_b  )
+                  CALL iom_get( numror, 'emp_corr'  , emp_corr )
+               ELSE                                                                                 !    as specified in namelist
+                  IF(lwp)   WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0'
+                  emp_corr = rn_fwb0
+                  a_fwb    = 999._wp
+                  a_fwb_b  = 999._wp    
+               END IF
                !
-               a_fwb_ini = a_fwb_b
-            ELSE                                                                                 !    as specified in namelist
-               IF(lwp)   WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0'
-               a_fwb   = rn_fwb0
-               a_fwb_b = 0._wp   ! used only the first year then it is replaced by a_fwb_ini
+               IF(lwp)   WRITE(numout,*)
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : initial freshwater correction flux = ', emp_corr , 'kg/m2/s'
                !
-               a_fwb_ini = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) &
-                  &      * rho0 / ( area * rday * REAL(nyear_len(1), wp) )
-            END IF
-            !
-            IF(lwp)   WRITE(numout,*)
-            IF(lwp)   WRITE(numout,*)'sbc_fwb : freshwater-budget at the end of previous year = ', a_fwb    , 'kg/m2/s'
-            IF(lwp)   WRITE(numout,*)'          freshwater-budget at initial state            = ', a_fwb_ini, 'kg/m2/s'
+            ENDIF
             !
-         ELSE
             ! at the end of year n:
             ikty = nyear_len(1) * 86400 / NINT(rn_Dt)
-            IF( MOD( kt, ikty ) == 0 ) THEN   ! Update a_fwb at the last time step of a year
-               !                                It should be the first time step of a year MOD(kt-1,ikty) but then the restart would be wrong
-               !                                Hence, we make a small error here but the code is restartable
-               a_fwb_b = a_fwb_ini
+            IF( MOD( kt-1, ikty ) == 0 ) THEN   ! Update a_fwb at the last time step of a year
+               a_fwb_b = a_fwb
                ! mean sea level taking into account ice+snow
-               a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
-               a_fwb   = a_fwb * rho0 / ( area * rday * REAL(nyear_len(1), wp) )   ! convert in kg/m2/s
+#if defined key_agrif
+               agrif_sum = 0.0_wp                           ! set work scalar to 0
+               CALL Agrif_step_child(glob_sum_volume_agrif) ! integrate over grids
+               a_fwb = agrif_sum
+#else                  
+               a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass_b(:,:) * r1_rho0 ) )
+#endif
+               a_fwb = a_fwb * rho0 / area
+               !
+               ! Special case if less than a year has been performed:
+               ! hence namelist rn_fwb0 still rules
+               IF ( a_fwb_b == 999._wp ) a_fwb_b = a_fwb
+               !
+               emp_corr = ( a_fwb - a_fwb_b ) / ( rday * REAL(nyear_len(0), wp) ) + emp_corr
+               IF(lwp)   WRITE(numout,*)
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : Compute new global mass at step = ', kt
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : New      averaged liquid height (ocean + snow + ice) = ',    a_fwb * r1_rho0, 'm'
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : Previous averaged liquid height (ocean + snow + ice) = ',  a_fwb_b * r1_rho0, 'm'
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : Implied freshwater-budget correction flux = ', emp_corr , 'kg/m2/s'
             ENDIF
-            !
+#if defined key_agrif
+         ELSE ! child grid if any
+            IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+                emp_corr  =  Agrif_parent(emp_corr) 
+            ENDIF
+#endif
          ENDIF
          !
-         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes using previous year budget minus initial state
-            zcoef = ( a_fwb - a_fwb_b )
-            emp(:,:) = emp(:,:) + zcoef * tmask(:,:,1)
-            qns(:,:) = qns(:,:) - zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes
+            emp(:,:) = emp(:,:) + emp_corr * tmask(:,:,1)
+            qns(:,:) = qns(:,:) - emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
             ! outputs
-            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -zcoef * rcp * sst_m(:,:) * tmask(:,:,1) )
-            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -zcoef * tmask(:,:,1) )
-         ENDIF
-         ! Output restart information
-         IF( lrst_oce ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
-            IF(lwp) WRITE(numout,*) '~~~~'
-            CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b )
-            CALL iom_rstput( kt, nitrst, numrow, 'a_fwb',   a_fwb   )
-         END IF
-         !
-         IF( kt == nitend .AND. lwp ) THEN
-            WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', a_fwb  , 'kg/m2/s'
-            WRITE(numout,*) '          freshwater-budget at initial state                    = ', a_fwb_b, 'kg/m2/s'
+            IF( iom_use('hflx_fwb_cea') )  CALL iom_put( 'hflx_fwb_cea', -emp_corr * rcp * sst_m(:,:) * tmask(:,:,1) )
+            IF( iom_use('vflx_fwb_cea') )  CALL iom_put( 'vflx_fwb_cea', -emp_corr * tmask(:,:,1) )
          ENDIF
+
+         IF ( Agrif_Root() ) THEN
+            ! Output restart information (root grid only)
+            IF( lrst_oce ) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
+               IF(lwp) WRITE(numout,*) '~~~~'
+               CALL iom_rstput( kt, nitrst, numrow, 'a_fwb',   a_fwb   )
+               CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b )
+               CALL iom_rstput( kt, nitrst, numrow, 'emp_corr',emp_corr)
+            END IF
+            !
+            IF( kt == nitend .AND. lwp ) THEN
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : Previous          year averaged liquid height (ocean + snow + ice) = ',    a_fwb * r1_rho0, 'm'
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : Previous previous year averaged liquid height (ocean + snow + ice) = ',  a_fwb_b * r1_rho0, 'm'
+               IF(lwp)   WRITE(numout,*)'sbc_fwb : freshwater-budget correction flux = ', emp_corr , 'kg/m2/s'
+            END IF
+         END IF 
          !
       CASE ( 3 )                             !==  global fwf set to zero and spread out over erp area  ==!
          !
@@ -252,5 +369,30 @@ CONTAINS
       !
    END SUBROUTINE sbc_fwb
 
+#if defined key_agrif
+   SUBROUTINE glob_sum_area_agrif()
+      !!---------------------------------------------------------------------
+      !!           ***  compute area with embedded zooms ***
+      !!----------------------------------------------------------------------
+
+      agrif_sum = agrif_sum + glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:))
+
+   END SUBROUTINE glob_sum_area_agrif   
+
+   SUBROUTINE glob_sum_volume_agrif()
+      !!---------------------------------------------------------------------
+      !!           ***  Compute volume with embedded zooms  ***
+      !!----------------------------------------------------------------------
+
+      IF ( nn_ice==2 ) THEN
+         agrif_sum = agrif_sum + glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ( ssh(:,:,Kmm_a) + snwice_mass_b(:,:) * r1_rho0 ))
+      ELSE
+         agrif_sum = agrif_sum + glob_sum( 'sbcfwb', e1e2t(:,:) * tmask_agrif(:,:) * ssh(:,:,Kmm_a) )
+      ENDIF
+
+   END SUBROUTINE glob_sum_volume_agrif 
+
+#endif
+
    !!======================================================================
 END MODULE sbcfwb

src/OCE/SBC/sbcice_cice.F90

@@ -175,7 +175,7 @@ CONTAINS
       ! there is no restart file.
       ! Values from a CICE restart file would overwrite this
       IF( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1.) 
+         CALL nemo2cice( ts(:,:,1,jp_tem,Kmm) , sst , 'T' , 1._wp) 
       ENDIF  
 #endif
 
@@ -207,10 +207,10 @@ CONTAINS
       fr_iu(:,:)=0.0
       fr_iv(:,:)=0.0
 
-      CALL cice2nemo(aice,fr_i, 'T', 1. )
+      CALL cice2nemo(aice,fr_i, 'T', 1._wp )
       IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
-            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
+            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1._wp )
          ENDDO
       ENDIF
 
@@ -224,8 +224,8 @@ CONTAINS
       CALL lbc_lnk( 'sbcice_cice', fr_iu , 'U', 1.0_wp,  fr_iv , 'V', 1.0_wp )
 
       ! set the snow+ice mass
-      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
-      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1._wp )
+      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1._wp )
       snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
       snwice_mass_b(:,:) = snwice_mass(:,:)
 
@@ -315,7 +315,7 @@ CONTAINS
             ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
                                  + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
          END_2D
-         CALL nemo2cice(ztmp,strax,'F', -1. )
+         CALL nemo2cice(ztmp,strax,'F', -1._wp )
 
 ! y comp of wind stress (CI_2)
 ! V point to F point
@@ -323,7 +323,7 @@ CONTAINS
             ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
                                  + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
          END_2D
-         CALL nemo2cice(ztmp,stray,'F', -1. )
+         CALL nemo2cice(ztmp,stray,'F', -1._wp )
 
 ! Surface downward latent heat flux (CI_5)
          IF(ksbc == jp_flx) THEN
@@ -349,7 +349,7 @@ CONTAINS
             END_2D
          ENDIF
          DO jl=1,ncat
-            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1._wp )
 
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
@@ -358,7 +358,7 @@ CONTAINS
             ELSE
                ztmp(:,:) = botmelt(:,:,jl)
             ENDIF
-            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1._wp )
 
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
@@ -367,7 +367,7 @@ CONTAINS
             ELSE
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
             ENDIF
-            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
+            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1._wp )
          ENDDO
 
       ELSE IF(ksbc == jp_blk) THEN
@@ -375,39 +375,39 @@ CONTAINS
 ! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself)
 ! x comp and y comp of atmosphere surface wind (CICE expects on T points)
          ztmp(:,:) = wndi_ice(:,:)
-         CALL nemo2cice(ztmp,uatm,'T', -1. )
+         CALL nemo2cice(ztmp,uatm,'T', -1._wp )
          ztmp(:,:) = wndj_ice(:,:)
-         CALL nemo2cice(ztmp,vatm,'T', -1. )
+         CALL nemo2cice(ztmp,vatm,'T', -1._wp )
          ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
-         CALL nemo2cice(ztmp,wind,'T', 1. )    ! Wind speed (m/s)
+         CALL nemo2cice(ztmp,wind,'T', 1._wp )    ! Wind speed (m/s)
          ztmp(:,:) = qsr_ice(:,:,1)
-         CALL nemo2cice(ztmp,fsw,'T', 1. )     ! Incoming short-wave (W/m^2)
+         CALL nemo2cice(ztmp,fsw,'T', 1._wp )     ! Incoming short-wave (W/m^2)
          ztmp(:,:) = qlw_ice(:,:,1)
-         CALL nemo2cice(ztmp,flw,'T', 1. )     ! Incoming long-wave (W/m^2)
+         CALL nemo2cice(ztmp,flw,'T', 1._wp )     ! Incoming long-wave (W/m^2)
          ztmp(:,:) = tatm_ice(:,:)
-         CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
-         CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
+         CALL nemo2cice(ztmp,Tair,'T', 1._wp )    ! Air temperature (K)
+         CALL nemo2cice(ztmp,potT,'T', 1._wp )    ! Potential temp (K)
 ! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows  
          ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )    
                                                ! Constant (101000.) atm pressure assumed
-         CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
+         CALL nemo2cice(ztmp,rhoa,'T', 1._wp )    ! Air density (kg/m^3)
          ztmp(:,:) = qatm_ice(:,:)
-         CALL nemo2cice(ztmp,Qa,'T', 1. )      ! Specific humidity (kg/kg)
+         CALL nemo2cice(ztmp,Qa,'T', 1._wp )      ! Specific humidity (kg/kg)
          ztmp(:,:)=10.0
-         CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
+         CALL nemo2cice(ztmp,zlvl,'T', 1._wp )    ! Atmos level height (m)
 
 ! May want to check all values are physically realistic (as in CICE routine 
 ! prepare_forcing)?
 
 ! Divide shortwave into spectral bands (as in prepare_forcing)
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
-         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
+         CALL nemo2cice(ztmp,swvdr,'T', 1._wp )             
          ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
-         CALL nemo2cice(ztmp,swvdf,'T', 1. )              
+         CALL nemo2cice(ztmp,swvdf,'T', 1._wp )              
          ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
-         CALL nemo2cice(ztmp,swidr,'T', 1. )
+         CALL nemo2cice(ztmp,swidr,'T', 1._wp )
          ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
-         CALL nemo2cice(ztmp,swidf,'T', 1. )
+         CALL nemo2cice(ztmp,swidf,'T', 1._wp )
 
       ENDIF
 
@@ -415,37 +415,37 @@ CONTAINS
 ! Ensure fsnow is positive (as in CICE routine prepare_forcing)
       IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit  
       ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)  
-      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
+      CALL nemo2cice(ztmp,fsnow,'T', 1._wp ) 
 
 ! Rainfall
       IF( iom_use('precip') )   CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
       ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
-      CALL nemo2cice(ztmp,frain,'T', 1. ) 
+      CALL nemo2cice(ztmp,frain,'T', 1._wp ) 
 
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
       nfrzmlt(:,:) = rho0 * rcp * e3t_m(:,:) * ( Tocnfrz-sst_m(:,:) ) / ( 2.0*dt )
 
       ztmp(:,:) = nfrzmlt(:,:)
-      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
+      CALL nemo2cice(ztmp,frzmlt,'T', 1._wp )
 
 ! SST  and SSS
 
-      CALL nemo2cice(sst_m,sst,'T', 1. )
-      CALL nemo2cice(sss_m,sss,'T', 1. )
+      CALL nemo2cice(sst_m,sst,'T', 1._wp )
+      CALL nemo2cice(sss_m,sss,'T', 1._wp )
 
 ! x comp and y comp of surface ocean current
 ! U point to F point
       DO_2D( 1, 1, 1, 0 )
          ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice(ztmp,uocn,'F', -1. )
+      CALL nemo2cice(ztmp,uocn,'F', -1._wp )
 
 ! V point to F point
       DO_2D( 1, 0, 1, 1 )
          ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice(ztmp,vocn,'F', -1. )
+      CALL nemo2cice(ztmp,vocn,'F', -1._wp )
 
       IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
           !
@@ -470,14 +470,14 @@ CONTAINS
          ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  )) * r1_e1u(ji,jj  )    &
             &               + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1)  ) * fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice( ztmp,ss_tltx,'F', -1. )
+      CALL nemo2cice( ztmp,ss_tltx,'F', -1._wp )
 
 ! T point to F point
       DO_2D( 1, 0, 1, 0 )
          ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj)) * r1_e2v(ji  ,jj)    &
             &               + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj)  ) *  fmask(ji,jj,1)
       END_2D
-      CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
+      CALL nemo2cice(ztmp,ss_tlty,'F', -1._wp )
       !
    END SUBROUTINE cice_sbc_in
 
@@ -499,7 +499,7 @@ CONTAINS
       ENDIF
       
 ! x comp of ocean-ice stress 
-      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
+      CALL cice2nemo(strocnx,ztmp1,'F', -1._wp )
       ss_iou(:,:)=0.0
 ! F point to U point
       DO_2D( 0, 0, 0, 0 )
@@ -508,7 +508,7 @@ CONTAINS
       CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1.0_wp )
 
 ! y comp of ocean-ice stress 
-      CALL cice2nemo(strocny,ztmp1,'F', -1. )
+      CALL cice2nemo(strocny,ztmp1,'F', -1._wp )
       ss_iov(:,:)=0.0
 ! F point to V point
 
@@ -527,8 +527,8 @@ CONTAINS
  
 ! Also need ice/ocean stress on T points so that taum can be updated 
 ! This interpolation is already done in CICE so best to use those values 
-      CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 
-      CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 
+      CALL cice2nemo(strocnxT,ztmp1,'T',-1._wp) 
+      CALL cice2nemo(strocnyT,ztmp2,'T',-1._wp) 
  
 ! Update taum with modulus of ice-ocean stress 
 ! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here 
@@ -551,11 +551,11 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       ENDIF
 
 #if defined key_cice4
-      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
-      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
+      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1._wp )
+      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1._wp )
 #else
-      CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
-      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
+      CALL cice2nemo(fresh_ai,ztmp1,'T', 1._wp )
+      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1._wp )
 #endif
 
 ! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
@@ -589,9 +589,9 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
 ! Now add in ice / snow related terms
 ! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
 #if defined key_cice4
-      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
+      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1._wp )
 #else
-      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
+      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1._wp )
 #endif
       qsr(:,:)=qsr(:,:)+ztmp1(:,:)
       CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1.0_wp )
@@ -601,9 +601,9 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       END_2D
 
 #if defined key_cice4
-      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
+      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1._wp )
 #else
-      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
+      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1._wp )
 #endif
       qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
 
@@ -611,10 +611,10 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
 
 ! Prepare for the following CICE time-step
 
-      CALL cice2nemo(aice,fr_i,'T', 1. )
+      CALL cice2nemo(aice,fr_i,'T', 1._wp )
       IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
-            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
+            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1._wp )
          ENDDO
       ENDIF
 
@@ -628,8 +628,8 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       CALL lbc_lnk( 'sbcice_cice', fr_iu , 'U', 1.0_wp, fr_iv , 'V', 1.0_wp )
 
       ! set the snow+ice mass
-      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
-      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1._wp )
+      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1._wp )
       snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
       snwice_mass_b(:,:) = snwice_mass(:,:)
       snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
@@ -661,16 +661,16 @@ taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2)
       !
       ! x and y comp of ice velocity
       !
-      CALL cice2nemo(uvel,u_ice,'F', -1. )
-      CALL cice2nemo(vvel,v_ice,'F', -1. )
+      CALL cice2nemo(uvel,u_ice,'F', -1._wp )
+      CALL cice2nemo(vvel,v_ice,'F', -1._wp )
       !
       ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out  
       !
       ! Snow and ice thicknesses (CO_2 and CO_3)
       !
       DO jl = 1, ncat
-         CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
-         CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
+         CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1._wp )
+         CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1._wp )
       END DO
       !
    END SUBROUTINE cice_sbc_hadgam

src/OCE/SBC/sbcmod.F90

@@ -75,6 +75,7 @@ MODULE sbcmod
    INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "single_precision_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
    !! $Id: sbcmod.F90 15372 2021-10-14 15:47:24Z davestorkey $
@@ -125,6 +126,10 @@ CONTAINS
 #endif
 #if ! defined key_si3
       IF( nn_ice == 2 )    nn_ice = 0  ! without key key_si3 you cannot use si3...
+#endif
+#if defined key_agrif
+      ! In case of an agrif zoom, the freshwater water budget is determined by parent:
+      IF (.NOT.Agrif_Root()) nn_fwb = Agrif_Parent(nn_fwb) 
 #endif
       !
       !
@@ -443,8 +448,8 @@ CONTAINS
             vtau(ji,jj) = vtau(ji,jj) * ( tauoc_wave(ji,jj) + tauoc_wave(ji,jj-1) ) * 0.5_wp
          END_2D
          !
-         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
-         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1._wp )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1._wp )
          !
          taum(:,:) = taum(:,:)*tauoc_wave(:,:)
          !
@@ -454,8 +459,8 @@ CONTAINS
       ELSEIF( ln_wave .AND. ln_taw ) THEN                  ! Wave stress reduction
          utau(:,:) = utau(:,:) - tawx(:,:) + twox(:,:)
          vtau(:,:) = vtau(:,:) - tawy(:,:) + twoy(:,:)
-         CALL lbc_lnk( 'sbcwave', utau, 'U', -1. )
-         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1. )
+         CALL lbc_lnk( 'sbcwave', utau, 'U', -1._wp )
+         CALL lbc_lnk( 'sbcwave', vtau, 'V', -1._wp )
          !
          DO_2D( 0, 0, 0, 0)
              taum(ji,jj) = sqrt((.5*(utau(ji-1,jj)+utau(ji,jj)))**2 + (.5*(vtau(ji,jj-1)+vtau(ji,jj)))**2)
@@ -465,7 +470,7 @@ CONTAINS
             &                                'If not requested select ln_taw=.false.' )
          !
       ENDIF
-      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1. )
+      CALL lbc_lnk( 'sbcmod', taum(:,:), 'T', 1._wp )
       !
       IF( ln_icebergs ) THEN  ! save pure stresses (with no ice-ocean stress) for use by icebergs
          utau_icb(:,:) = utau(:,:) ; vtau_icb(:,:) = vtau(:,:) 
@@ -603,16 +608,16 @@ CONTAINS
       ENDIF
       !
       IF(sn_cfctl%l_prtctl) THEN     ! print mean trends (used for debugging)
-         CALL prt_ctl(tab2d_1=fr_i                , clinfo1=' fr_i     - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=(emp-rnf)           , clinfo1=' emp-rnf  - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=(sfx-rnf)           , clinfo1=' sfx-rnf  - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=qns                 , clinfo1=' qns      - : ', mask1=tmask )
-         CALL prt_ctl(tab2d_1=qsr                 , clinfo1=' qsr      - : ', mask1=tmask )
-         CALL prt_ctl(tab3d_1=tmask               , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
+         CALL prt_ctl(tab2d_1=CASTDP(fr_i)                , clinfo1=' fr_i     - : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP((emp-rnf))           , clinfo1=' emp-rnf  - : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP((sfx-rnf))           , clinfo1=' sfx-rnf  - : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(qns)                 , clinfo1=' qns      - : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=CASTDP(qsr)                , clinfo1=' qsr      - : ', mask1=tmask )
+         CALL prt_ctl(tab3d_1=CASTDP(tmask)              , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
          CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst      - : ', mask1=tmask, kdim=1   )
          CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss      - : ', mask1=tmask, kdim=1   )
-         CALL prt_ctl(tab2d_1=utau                , clinfo1=' utau     - : ', mask1=umask,                      &
-            &         tab2d_2=vtau                , clinfo2=' vtau     - : ', mask2=vmask )
+         CALL prt_ctl(tab2d_1=CASTDP(utau)               , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=CASTDP(vtau)               , clinfo2=' vtau     - : ', mask2=vmask )
       ENDIF
 
       IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary

src/OCE/SBC/sbcrnf.F90

@@ -20,6 +20,7 @@ MODULE sbcrnf
    USE sbc_oce        ! surface boundary condition variables
    USE eosbn2         ! Equation Of State
    USE closea, ONLY: l_clo_rnf, clo_rnf ! closed seas
+   USE isf_oce
    !
    USE in_out_manager ! I/O manager
    USE fldread        ! read input field at current time step
@@ -40,8 +41,10 @@ MODULE sbcrnf
    LOGICAL                    ::      ln_rnf_depth_ini  !: depth       river runoffs  computed at the initialisation
    REAL(wp)                   ::      rn_rnf_max        !: maximum value of the runoff climatologie (ln_rnf_depth_ini =T)
    REAL(wp)                   ::      rn_dep_max        !: depth over which runoffs is spread       (ln_rnf_depth_ini =T)
+   REAL(wp)                   ::      tot_flux           !: total iceberg flux (temporary variable)
    INTEGER                    ::      nn_rnf_depth_file !: create (=1) a runoff depth file or not (=0)
    LOGICAL           , PUBLIC ::   ln_rnf_icb        !: iceberg flux is specified in a file
+   LOGICAL                    ::   ln_icb_mass       !: Scale iceberg flux to match FW flux from coupled model
    LOGICAL                    ::   ln_rnf_tem        !: temperature river runoffs attribute specified in a file
    LOGICAL           , PUBLIC ::   ln_rnf_sal        !: salinity    river runoffs attribute specified in a file
    TYPE(FLD_N)       , PUBLIC ::   sn_rnf            !: information about the runoff file to be read
@@ -119,8 +122,8 @@ CONTAINS
       !
       IF( .NOT. l_rnfcpl )  THEN
                             CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt ( runoffs + iceberg )
-         IF( ln_rnf_icb )   CALL fld_read ( kt, nn_fsbc, sf_i_rnf )    ! idem for iceberg flux if required
       ENDIF
+      IF(   ln_rnf_icb   )   CALL fld_read ( kt, nn_fsbc, sf_i_rnf )    ! idem for iceberg flux        if required
       IF(   ln_rnf_tem   )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required
       IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
       !
@@ -138,14 +141,57 @@ CONTAINS
                 CALL iom_put( 'hflx_icb_cea' , -fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
              ENDIF
          ENDIF
-         !
+
+         IF( ln_rnf_icb ) THEN
+            fwficb(:,:) = rn_rfact * ( sf_i_rnf(1)%fnow(:,:,1) ) * tmask(:,:,1)  ! updated runoff value at time step kt
+            IF( ln_icb_mass ) THEN
+                ! Modify the Iceberg FW flux to be consistent with the change in
+                ! mass of the Antarctic/Greenland ice sheet for FW conservation in
+                ! coupled model. This isn't perfect as FW flux will go into ocean at
+                ! wrong time of year but more important to maintain FW balance
+                tot_flux = SUM(fwficb(:,:)*e1e2t(:,:)*tmask_i(:,:)*greenland_icesheet_mask(:,:)) !Need to multiply by area to convert to kg/s
+                IF( lk_mpp ) CALL mpp_sum( 'icbclv', tot_flux )
+                IF( tot_flux > rsmall ) THEN
+                    WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                                 &
+                    &    fwficb(:,:) = fwficb(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction &
+                    &                                     / tot_flux
+	                ELSE IF( rn_greenland_calving_fraction < rsmall ) THEN
+                    WHERE( greenland_icesheet_mask(:,:) == 1.0 ) fwficb(:,:) = 0.0
+                ELSE
+                    CALL CTL_STOP('STOP', 'No iceberg runoff data read in for Greenland. Check input file or set rn_greenland_calving_fraction=0.0')
+                ENDIF
+
+                tot_flux = SUM(fwficb(:,:)*e1e2t(:,:)*tmask_i(:,:)*antarctica_icesheet_mask(:,:))
+                IF( lk_mpp ) CALL mpp_sum( 'icbclv', tot_flux )
+                IF( tot_flux > rsmall ) THEN
+                    WHERE( antarctica_icesheet_mask(:,:) == 1.0 )                                                                                &
+                    &    fwficb(:,:) = fwficb(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction &
+                    &                                     / (tot_flux  + 1.0e-10_wp )
+                ELSE IF( rn_antarctica_calving_fraction < rsmall ) THEN
+                    WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) fwficb(:,:) = 0.0
+                ELSE
+                    CALL CTL_STOP('STOP', 'No iceberg runoff data read in for Antarctica. Check input file or set rn_antarctica_calving_fraction=0.0')
+                ENDIF
+            ENDIF
+            CALL lbc_lnk('sbcrnf',rnf,'T',1._wp,fwficb,'T',1._wp)
+            CALL iom_put( 'berg_melt'  , fwficb(:,:)  )          ! output iceberg flux
+            CALL iom_put( 'hflx_icb_cea' , fwficb(:,:) * rLfus )   ! output Heat Flux into Sea Water due to Iceberg Thermodynamics -->
+            rnf(:,:) = rnf(:,:) + fwficb(:,:)                      ! fwficb isn't used anywhere else so add it to runoff here
+            qns_tot(:,:) = qns_tot(:,:) - fwficb(:,:) * rLfus      ! TG: I think this is correct
+         ENDIF         !
+
          !                                                           ! set temperature & salinity content of runoffs
          IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
             rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rho0
             CALL eos_fzp( sss_m(:,:), ztfrz(:,:) )
             WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
-               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rho0
+               rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
             END WHERE
+            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )
+            ! where fwf comes from melting of ice shelves or iceberg
+                rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rLfusisf * r1_rho0_rcp
+            END WHERE
+
          ELSE                                                        ! use SST as runoffs temperature
             !CEOD River is fresh water so must at least be 0 unless we consider ice
             rnf_tsc(:,:,jp_tem) = MAX( sst_m(:,:), 0.0_wp ) * rnf(:,:) * r1_rho0
@@ -158,11 +204,15 @@ CONTAINS
          IF( iom_use('sflx_rnf_cea') )   CALL iom_put( 'sflx_rnf_cea', rnf_tsc(:,:,jp_sal) * rho0       )   ! output runoff salt flux (g/m2/s)
       ENDIF
       !
+
+      ! Make sure rnf is up to date!
+      call lbc_lnk('rnf_div', rnf, 'T', 1.0)
+
       !                                                ! ---------------------------------------- !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
          !                                             ! ---------------------------------------- !
          IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN         !* Restart: read in restart file
-            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields red in the restart file', lrxios
+            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields read in the restart file', lrxios
             CALL iom_get( numror, jpdom_auto, 'rnf_b'   , rnf_b                 )   ! before runoff
             CALL iom_get( numror, jpdom_auto, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem) )   ! before heat content of runoff
             CALL iom_get( numror, jpdom_auto, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal) )   ! before salinity content of runoff
@@ -259,7 +309,7 @@ CONTAINS
       REAL(wp)          ::   zacoef
       REAL(wp), DIMENSION(jpi,jpj,2) :: zrnfcl
       !!
-      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb,   &
+      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal, ln_rnf_icb, ln_icb_mass, &
          &                 sn_rnf, sn_cnf    , sn_i_rnf, sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
          &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact,     &
          &                 ln_rnf_depth_ini  , rn_dep_max  , rn_rnf_max, nn_rnf_depth_file
@@ -268,16 +318,6 @@ CONTAINS
       !                                         !==  allocate runoff arrays
       IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' )
       !
-      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths
-         ln_rnf_mouth  = .FALSE.                   ! default definition needed for example by sbc_ssr or by tra_adv_muscl
-         nkrnf         = 0
-         rnf     (:,:) = 0.0_wp
-         rnf_b   (:,:) = 0.0_wp
-         rnfmsk  (:,:) = 0.0_wp
-         rnfmsk_z(:)   = 0.0_wp
-         RETURN
-      ENDIF
-      !
       !                                   ! ============
       !                                   !   Namelist
       !                                   ! ============
@@ -289,6 +329,19 @@ CONTAINS
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_rnf in configuration namelist' )
       IF(lwm) WRITE ( numond, namsbc_rnf )
       !
+      IF( .NOT. ln_rnf ) THEN                      ! no specific treatment in vicinity of river mouths
+         ln_rnf_mouth  = .FALSE.                   ! default definition needed for example by sbc_ssr or by tra_adv_muscl
+         ln_rnf_tem    = .FALSE.
+         ln_rnf_sal    = .FALSE.
+         ln_rnf_icb    = .FALSE.
+         nkrnf         = 0
+         rnf     (:,:) = 0.0_wp
+         rnf_b   (:,:) = 0.0_wp
+         rnfmsk  (:,:) = 0.0_wp
+         rnfmsk_z(:)   = 0.0_wp
+         RETURN
+      ENDIF
+      !
       !                                         ! Control print
       IF(lwp) THEN
          WRITE(numout,*)
@@ -315,21 +368,22 @@ CONTAINS
          IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
          CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf', no_print )
          !
-         IF( ln_rnf_icb ) THEN                      ! Create (if required) sf_i_rnf structure
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '          iceberg flux read in a file'
-            ALLOCATE( sf_i_rnf(1), STAT=ierror  )
-            IF( ierror > 0 ) THEN
-               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' )   ;   RETURN
-            ENDIF
-            ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1)   )
-            IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
-         ELSE
-            fwficb(:,:) = 0._wp
-         ENDIF
 
       ENDIF
+
+      IF( ln_rnf_icb ) THEN                      ! Create (if required) sf_i_rnf structure
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          iceberg flux read in a file'
+         ALLOCATE( sf_i_rnf(1), STAT=ierror  )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_i_rnf structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_i_rnf(1)%fnow(jpi,jpj,1)   )
+         IF( sn_i_rnf%ln_tint ) ALLOCATE( sf_i_rnf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill (sf_i_rnf, (/ sn_i_rnf /), cn_dir, 'sbc_rnf_init', 'read iceberg flux data', 'namsbc_rnf' )
+      ELSE
+         fwficb(:,:) = 0._wp
+      ENDIF
       !
       IF( ln_rnf_tem ) THEN                      ! Create (if required) sf_t_rnf structure
          IF(lwp) WRITE(numout,*)

src/OCE/SBC/sbcssm.F90

@@ -31,6 +31,7 @@ MODULE sbcssm
 
    LOGICAL, SAVE ::   l_ssm_mean = .FALSE.   ! keep track of whether means have been read from restart file
 
+#  include "single_precision_substitute.h90"
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -58,13 +59,28 @@ CONTAINS
       INTEGER  ::   ji, jj               ! loop index
       REAL(wp) ::   zcoef, zf_sbc       ! local scalar
       REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts
-      !!---------------------------------------------------------------------
+      CHARACTER(len=4),SAVE :: stype 
+      !!--------------------------------------------------------------------- 
+      IF( kt == nit000 ) THEN 
+         IF( ln_TEOS10 ) THEN 
+            stype='abs'    ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module)
+         ELSE IF( ln_SEOS ) THEN
+            stype='seos'   ! seos: using Simplified Equation of state (sst is converted to potential temperature for the surface module)
+         ELSE
+            stype='pra'    ! eos-80: using practical salinity
+         ENDIF
+      ENDIF 
       !
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
       zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm)
       zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
       !
-         !                                                ! ---------------------------------------- !
+      !                                       !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
+      !                                       !                otherwise restartability and reproducibility are broken 
+      !                                       !                computed in traqsr only on the inner domain 
+      CALL lbc_lnk( 'sbc_ssm', fraqsr_1lev(:,:), 'T', 1._wp )
+      !
+      !                                                   ! ---------------------------------------- !
       IF( nn_fsbc == 1 ) THEN                             !      Instantaneous surface fields        !
          !                                                ! ---------------------------------------- !
          ssu_m(:,:) = uu(:,:,1,Kbb)
@@ -170,8 +186,8 @@ CONTAINS
       IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   !
          CALL iom_put( 'ssu_m', ssu_m )
          CALL iom_put( 'ssv_m', ssv_m )
-         CALL iom_put( 'sst_m', sst_m )
-         CALL iom_put( 'sss_m', sss_m )
+         CALL iom_put( 'sst_m_pot', sst_m ) 
+         CALL iom_put( 'sss_m_'//stype, sss_m ) 
          CALL iom_put( 'ssh_m', ssh_m )
          CALL iom_put( 'e3t_m', e3t_m )
          CALL iom_put( 'frq_m', frq_m )
@@ -241,7 +257,7 @@ CONTAINS
          IF(lwp) WRITE(numout,*) '   default initialisation of ss._m arrays'
          ssu_m(:,:) = uu(:,:,1,Kbb)
          ssv_m(:,:) = vv(:,:,1,Kbb)
-         IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( ts(:,:,1,jp_tem,Kmm), ts(:,:,1,jp_sal,Kmm) )
+         IF( l_useCT )  THEN    ;   sst_m(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
          ELSE                   ;   sst_m(:,:) = ts(:,:,1,jp_tem,Kmm)
          ENDIF
          sss_m(:,:) = ts  (:,:,1,jp_sal,Kmm)

src/OCE/SBC/sbcssr.F90

@@ -93,15 +93,14 @@ CONTAINS
          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN      !    Add restoring term     !
             !                                      ! ========================= !
             !
-            qrp(:,:) = 0._wp ! necessary init
-            erp(:,:) = 0._wp
-            !
-            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term
+            IF(     nn_sstr == 1 ) THEN                                   !* Temperature restoring term
                DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
                   zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1)
                   qns(ji,jj) = qns(ji,jj) + zqrp
                   qrp(ji,jj) = zqrp
                END_2D
+            ELSEIF( nn_sssr == 2 ) THEN
+               qrp(:,:) = 0._wp   ! necessary init, see bellow: qrp(ji,jj) = qrp(ji,jj) - ...
             ENDIF
             !
             IF( nn_sssr /= 0 .AND. nn_sssr_ice /= 1 ) THEN
@@ -233,8 +232,8 @@ CONTAINS
       !
       coefice(:,:) = 1._wp         !  Initialise coefice to 1._wp ; will not need to be changed if nn_sssr_ice=1
       !                            !* Initialize qrp and erp if no restoring 
-      IF( nn_sstr /= 1                   )   qrp(:,:) = 0._wp
-      IF( nn_sssr /= 1 .OR. nn_sssr /= 2 )   erp(:,:) = 0._wp
+      IF( nn_sstr /= 1 .AND. nn_sssr /= 2 )   qrp(:,:) = 0._wp
+      IF( nn_sssr /= 1 .AND. nn_sssr /= 2 )   erp(:,:) = 0._wp
       !
    END SUBROUTINE sbc_ssr_init
          

src/OCE/SBC/sbcwave.F90

@@ -71,7 +71,8 @@ MODULE sbcwave
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tsd2d           !: Surface Stokes Drift module at t-point
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   div_sd          !: barotropic stokes drift divergence
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   ut0sd, vt0sd    !: surface Stokes drift velocities at t-point
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   usd, vsd, wsd   !: Stokes drift velocities at u-, v- & w-points, resp.u
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  :: wsd!: Stokes drift velocities at u-, v- & w-points, resp.u
+   REAL(dp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  :: usd, vsd!: Stokes drift velocities at u-, v- & w-points, resp.u
 !
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   charn           !: charnock coefficient at t-point
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   tawx            !: Net wave-supported stress, u
@@ -211,7 +212,7 @@ CONTAINS
          END_3D
       ENDIF
 
-      CALL lbc_lnk( 'sbcwave', usd, 'U', -1.0_wp, vsd, 'V', -1.0_wp )
+      CALL lbc_lnk( 'sbcwave', usd, 'U', -1.0_dp, vsd, 'V', -1.0_dp )
 
       !
       !                       !==  vertical Stokes Drift 3D velocity  ==!