Changeset 9987 for branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC

Timestamp:

2018-07-23T11:33:03+02:00 (6 years ago)

Author:

emmafiedler

Message:

Merge with GO6 FOAMv14 package branch r9288

Location:

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC

Files:

• albedo.F90 (modified) (7 diffs)
• cpl_oasis3.F90 (modified) (11 diffs)
• fldread.F90 (modified) (4 diffs)
• geo2ocean.F90 (modified) (2 diffs)
• sbc_ice.F90 (modified) (4 diffs)
• sbc_oce.F90 (modified) (3 diffs)
• sbcblk_clio.F90 (modified) (1 diff)
• sbcblk_core.F90 (modified) (7 diffs)
• sbccpl.F90 (modified) (54 diffs)
• sbcfwb.F90 (modified) (2 diffs)
• sbcice_cice.F90 (modified) (18 diffs)
• sbcice_if.F90 (modified) (1 diff)
• sbcice_lim.F90 (modified) (10 diffs)
• sbcice_lim_2.F90 (modified) (1 diff)
• sbcisf.F90 (modified) (16 diffs)
• sbcmod.F90 (modified) (4 diffs)
• sbcrnf.F90 (modified) (3 diffs)
• updtide.F90 (modified) (5 diffs)

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

@@ -9 +9 @@
    !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0
    !!             -   ! 2006-08  (G. Madec)  cleaning for surface module
+   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo
    !!----------------------------------------------------------------------
 
@@ -29 +30 @@
 
    INTEGER  ::   albd_init = 0      !: control flag for initialization
-   REAL(wp) ::   zzero     = 0.e0   ! constant values
-   REAL(wp) ::   zone      = 1.e0   !    "       "
-
-   REAL(wp) ::   c1     = 0.05    ! constants values
-   REAL(wp) ::   c2     = 0.10    !    "        "
-   REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude
-
+  
+   REAL(wp) ::   rmue     = 0.40    !  cosine of local solar altitude
+   REAL(wp) ::   ralb_oce = 0.066   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   REAL(wp) ::   c1       = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) ::   c2       = 0.10    !  "        "
+   REAL(wp) ::   rcloud   = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+ 
    !                             !!* namelist namsbc_alb
-   REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984)
-#if defined key_lim3
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#else
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#endif
-   REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute
-   REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972)
-   REAL(wp) ::   rn_alphc         ! 
+   INTEGER  ::   nn_ice_alb
+   REAL(wp) ::   rn_albice
 
    !!----------------------------------------------------------------------
@@ -59 +53 @@
       !!          
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!                as well as the ocean one
       !!       
-      !! ** Method  : - Computation of the albedo of snow or ice (choose the 
-      !!                rignt one by a large number of tests
-      !!              - Computation of the albedo of the ocean
-      !!
-      !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
+      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
+      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                           and Grenfell & Perovich (JGR 2004)
+      !!                Description of scheme 1:
+      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
+      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
+      !!                     0-5cm  : linear function of ice thickness
+      !!                     5-150cm: log    function of ice thickness
+      !!                     > 150cm: constant
+      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
+      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
+      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
+      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
+      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
+      !!
+      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                     under melting conditions than under freezing conditions
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
+      !!
+      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !!                Brandt et al. 2005, J. Climate, vol 18
+      !!                Grenfell & Perovich 2004, JGR, vol 109 
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -73 +87 @@
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
       !!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays)
-      REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting
-      REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing
-      REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow
-      REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow
-      REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
-      REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow)
-      REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1
-      REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2
-      !!
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      INTEGER  ::   ijpl               ! number of ice categories (3rd dim of ice input arrays)
+      REAL(wp)            ::   ralb_im, ralb_sf, ralb_sm, ralb_if
+      REAL(wp)            ::   zswitch, z1_c1, z1_c2
+      REAL(wp)                            ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it             ! intermediate variable & albedo of ice (snow free)
       !!---------------------------------------------------------------------
-      
+
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-
-      CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth )
+      
+      CALL wrk_alloc( jpi,jpj,ijpl, zalb, zalb_it )
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
 
-      !---------------------------
-      !  Computation of  zficeth
-      !---------------------------
-      ! ice free of snow and melts
-      WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice
-      ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
-      END  WHERE
-
-      WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz
-      ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
-      ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
-         &                                                                 - 0.8608 * ph_ice * ph_ice     &
-         &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-      ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
-      END WHERE
-
-!!gm old code
-!      DO jl = 1, ijpl
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl)
-!               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 )
-!               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN
-!                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               &
-!                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 &
-!                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl)
-!               ELSE
-!                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl) 
-!               ENDIF
-!            END DO
-!         END DO
-!      END DO
-!!gm end old code
-      
-      !----------------------------------------------- 
-      !    Computation of the snow/ice albedo system 
-      !-------------------------- ---------------------
-      
-      !    Albedo of snow-ice for clear sky.
-      !-----------------------------------------------    
-      DO jl = 1, ijpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !  Case of ice covered by snow.             
-               !                                        !  freezing snow        
-               zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) )
-               zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             &
-                  &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   &
-                  &     +         zihsc1   * rn_alphd  
-               !                                        !  melting snow                
-               zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) )
-               zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   &
-                  &     +         zihsc2   *   rn_alphc 
-               !
-               zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )   
-               zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf
-            
-               !  Case of ice free of snow.
-               zalbpic  = zficeth(ji,jj,jl) 
-            
-               ! albedo of the system   
-               zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) )
-               pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic
+      
+      SELECT CASE ( nn_ice_alb )
+
+      !------------------------------------------
+      !  Shine and Henderson-Sellers (1985)
+      !------------------------------------------
+      CASE( 0 )
+       
+         ralb_sf = 0.80       ! dry snow
+         ralb_sm = 0.65       ! melting snow
+         ralb_if = 0.72       ! bare frozen ice
+         ralb_im = rn_albice  ! bare puddled ice 
+         
+         !  Computation of ice albedo (free of snow)
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
+         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
+         END  WHERE
+      
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
+            &                                                                 - 0.8608 * ph_ice * ph_ice     &
+            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
+         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
+         END WHERE
+     
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! freezing snow
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
+                  !                                        !  freezing snow        
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
+                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
+                     &        +         zswitch   * ralb_sf  
+
+                  ! melting snow
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
+                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
+                      &     +         zswitch   *   ralb_sm 
+                  !
+                  ! snow albedo
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+               
+                  ! Ice/snow albedo
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
+                  !
+               END DO
             END DO
          END DO
-      END DO
-      
-      !    Albedo of snow-ice for overcast sky.
-      !----------------------------------------------  
-      pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction
-      !
-      CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth )
+
+         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+
+      !------------------------------------------
+      !  New parameterization (2016)
+      !------------------------------------------
+      CASE( 1 ) 
+
+         ralb_im = rn_albice  ! bare puddled ice
+! compilation of values from literature
+         ralb_sf = 0.85      ! dry snow
+         ralb_sm = 0.75      ! melting snow
+         ralb_if = 0.60      ! bare frozen ice
+! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
+!         ralb_sf = 0.85       ! dry snow
+!         ralb_sm = 0.72       ! melting snow
+!         ralb_if = 0.65       ! bare frozen ice
+! Brandt et al 2005 (East Antarctica)
+!         ralb_sf = 0.87      ! dry snow
+!         ralb_sm = 0.82      ! melting snow
+!         ralb_if = 0.54      ! bare frozen ice
+! 
+         !  Computation of ice albedo (free of snow)
+         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) ) 
+         z1_c2 = 1. / 0.05
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
+         ELSE WHERE                                              ;   zalb = ralb_if
+         END  WHERE
+         
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18 - zalb     ) * z1_c1 *  &
+            &                                                                     ( LOG(1.5) - LOG(ph_ice) )
+         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18 - ralb_oce ) * z1_c2 * ph_ice
+         END WHERE
+
+         z1_c1 = 1. / 0.02
+         z1_c2 = 1. / 0.03
+         !  Computation of the snow/ice albedo
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
+                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
+
+                   ! snow albedo
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+
+                  ! Ice/snow albedo   
+                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl)
+
+              END DO
+            END DO
+         END DO
+         ! Effect of the clouds (2d order polynomial)
+         pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ); 
+
+      END SELECT
+      
+      CALL wrk_dealloc( jpi,jpj,ijpl, zalb, zalb_it )
       !
    END SUBROUTINE albedo_ice
@@ -181 +229 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) ::   zcoef   ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982 
-      pa_oce_cs(:,:) = zcoef               
-      pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      REAL(wp) :: zcoef 
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
+      pa_oce_cs(:,:) = zcoef 
+      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
    END SUBROUTINE albedo_oce
@@ -200 +248 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
       !!----------------------------------------------------------------------
       !
@@ -219 +267 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
-         WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud
-         WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice
-         WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd
-         WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi
-         WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc
+         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
+         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
       ENDIF
       !

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90

@@ -31 +31 @@
    USE in_out_manager               ! I/O manager
    USE lbclnk                       ! ocean lateral boundary conditions (or mpp link)
-
+   
    IMPLICIT NONE
    PRIVATE
@@ -41 +41 @@
    PUBLIC   cpl_freq
    PUBLIC   cpl_finalize
+#if defined key_mpp_mpi
+   INCLUDE 'mpif.h'
+#endif
+   
+   INTEGER, PARAMETER         :: localRoot  = 0
+   LOGICAL                    :: commRank            ! true for ranks doing OASIS communication
+#if defined key_cpl_rootexchg
+   LOGICAL                    :: rootexchg =.true.   ! logical switch 
+#else
+   LOGICAL                    :: rootexchg =.false.  ! logical switch 
+#endif 
 
    INTEGER, PUBLIC            ::   OASIS_Rcv  = 1    !: return code if received field
@@ -82 +93 @@
 
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   exfld   ! Temporary buffer for receiving
-
+   INTEGER, PUBLIC :: localComm 
+      
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -120 +132 @@
       IF ( nerror /= OASIS_Ok ) &
          CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' )
+      localComm = kl_comm 
       !
    END SUBROUTINE cpl_init
@@ -177 +190 @@
       IF( nerror > 0 ) THEN
          CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN
-      ENDIF
+      ENDIF      
       !
       ! -----------------------------------------------------------------
       ! ... Define the partition 
       ! -----------------------------------------------------------------
-      
+            
       paral(1) = 2                                              ! box partitioning
       paral(2) = jpiglo * (nldj-1+njmpp-1) + (nldi-1+nimpp-1)   ! NEMO lower left corner global offset    
@@ -196 +209 @@
       ENDIF
       
-      CALL oasis_def_partition ( id_part, paral, nerror )
+      CALL oasis_def_partition ( id_part, paral, nerror, jpiglo*jpjglo)
       !
       ! ... Announce send variables. 
@@ -241 +254 @@
             END DO
          ENDIF
-      END DO
+      END DO      
       !
       ! ... Announce received variables. 
@@ -373 +386 @@
             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, kinfo )   
                
                llaction =  kinfo == OASIS_Recvd   .OR. kinfo == OASIS_FromRest .OR.   &
@@ -384 +397 @@
                   kinfo = OASIS_Rcv
                   IF( llfisrt ) THEN 
-                     pdata(nldi:nlei,nldj:nlej,jc) =                                 exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm)
+                     pdata(nldi:nlei,nldj:nlej,jc) =                                 exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm) 
                      llfisrt = .FALSE.
                   ELSE
@@ -463 +476 @@
          CALL oasis_get_freqs(id, mop, 1, itmp, info)
 #else
+#if defined key_oasis3 
+         itmp(1) = namflddti( id )
+#else
          CALL oasis_get_freqs(id,      1, itmp, info)
+#endif
 #endif
          cpl_freq = itmp(1)
@@ -514 +531 @@
    END SUBROUTINE oasis_get_localcomm
 
-   SUBROUTINE oasis_def_partition(k1,k2,k3)
+   SUBROUTINE oasis_def_partition(k1,k2,k3,K4)
       INTEGER     , INTENT(  out) ::  k1,k3
       INTEGER     , INTENT(in   ) ::  k2(5)
+      INTEGER     , OPTIONAL, INTENT(in   ) ::  k4
       k1 = k2(1) ; k3 = k2(5)
       WRITE(numout,*) 'oasis_def_partition: Error you sould not be there...'

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

@@ -32 +32 @@
    PUBLIC   fld_map    ! routine called by tides_init
    PUBLIC   fld_read, fld_fill   ! called by sbc... modules
+   PUBLIC   fld_clopn
 
    TYPE, PUBLIC ::   FLD_N      !: Namelist field informations
@@ -815 +816 @@
          imonth = kmonth
          iday = kday
+         IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
+            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 )  
+            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month
+            llprevyr   = llprevmth .AND. nmonth == 1
+            iyear  = nyear  - COUNT((/llprevyr /))
+            imonth = nmonth - COUNT((/llprevmth/)) + 12 * COUNT((/llprevyr /))
+            iday   = nday   + nmonth_len(nmonth-1) * COUNT((/llprevmth/)) - isec_week / NINT(rday)
+         ENDIF
       ELSE                                                  ! use current day values
          IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week
@@ -1281 +1290 @@
       CHARACTER(LEN=*)          , INTENT(in   ) ::   lsmfile ! land sea mask file name
       !! 
-      REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta,zfieldo                  ! temporary array of values on input grid
+      REAL(wp),DIMENSION(:,:,:),ALLOCATABLE     ::   ztmp_fly_dta                          ! temporary array of values on input grid
       INTEGER, DIMENSION(3)                     ::   rec1,recn                             ! temporary arrays for start and length
       INTEGER, DIMENSION(3)                     ::   rec1_lsm,recn_lsm                     ! temporary arrays for start and length in case of seaoverland
@@ -1347 +1356 @@
 
 
-         itmpi=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),1)
-         itmpj=SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),2)
+         itmpi=jpi2_lsm-jpi1_lsm+1
+         itmpj=jpj2_lsm-jpj1_lsm+1
          itmpz=kk
          ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90

@@ -51 +51 @@
 
    SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
-                       px2 , py2 )
+                       px2 , py2 , kchoix  )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE repcmo  ***
@@ -68 +68 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point)
       !!----------------------------------------------------------------------
-      
-      ! Change from geographic to stretched coordinate
-      ! ----------------------------------------------
-      CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
-      CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
-      
+      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
 

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -80 +80 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_oce       !: heat flux of precip and evap over ocean     [W/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qemp_ice       !: heat flux of precip and evap over ice       [W/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: heat flux of precip over ice                [J/m3]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qevap_ice      !: heat flux of evap over ice                  [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qprec_ice      !: enthalpy of precip over ice                 [J/m3]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s]
 #endif
@@ -101 +102 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point
-   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz             !: sea surface freezing temperature
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice           !: sea-ice surface skin temperature (on categories)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   kn_ice             !: sea-ice surface layer thermal conductivity (on cats)
+
    ! variables used in the coupled interface
    INTEGER , PUBLIC, PARAMETER ::   jpl = ncat
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_p, ht_p ! Meltpond fraction and depth
+   
+   !
+   
+   !
+#if defined key_asminc
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ndaice_da          !: NEMO fresh water flux to ocean due to data assim
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfresh_da          !: NEMO salt flux to ocean due to data assim
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nfsalt_da          !: NEMO ice concentration change/second from data assim
+#endif
+      
 #endif
    
@@ -144 +159 @@
 #endif
 #if defined key_lim3
-         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,  &
-         &      qemp_ice(jpi,jpj)     , qemp_oce(jpi,jpj)      ,                       &
-         &      qns_oce (jpi,jpj)     , qsr_oce (jpi,jpj)      , emp_oce (jpi,jpj)  ,  &
+         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   &
+         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   &
+         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce (jpi,jpj)  ,   &
 #endif
          &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) )
@@ -152 +167 @@
 
 #if defined key_cice
-      ALLOCATE( qla_ice(jpi,jpj,1)    , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , &
+      ALLOCATE( qla_ice(jpi,jpj,ncat) , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , &
                 wndi_ice(jpi,jpj)     , tatm_ice(jpi,jpj)     , qatm_ice(jpi,jpj)     , &
                 wndj_ice(jpi,jpj)     , nfrzmlt(jpi,jpj)      , ss_iou(jpi,jpj)       , &
                 ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , &
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
-                STAT= ierr(1) )
-      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , &
+#if defined key_asminc
+                ndaice_da(jpi,jpj)    , nfresh_da(jpi,jpj)    , nfsalt_da(jpi,jpj)    , &
+#endif
+                sstfrz(jpi,jpj)       , STAT= ierr(1) )
+   ! Alex West: Allocating tn_ice with 5 categories.  When NEMO is used with CICE, this variable
+   ! represents top layer ice temperature, which is multi-category.
+      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,jpl)  , &
          &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , &
          &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , &
-         &                     STAT= ierr(2) )
+         &                     a_p(jpi,jpj,jpl)      , ht_p(jpi,jpj,jpl)     , tsfc_ice(jpi,jpj,jpl) , &
+         &                     kn_ice(jpi,jpj,jpl) ,    STAT=ierr(2) )
       
 #endif

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

@@ -125 +125 @@
 #endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_icesheet_mass_array, greenland_icesheet_mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mass_array, antarctica_icesheet_mask
 
    !!----------------------------------------------------------------------
@@ -137 +139 @@
    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. 
 
    !! * Substitutions
@@ -172 +201 @@
          &      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_mass_array(jpi,jpj) , antarctica_icesheet_mass_array(jpi,jpj) )
+      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) )
          !
 #if defined key_vvl

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

@@ -684 +684 @@
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
 
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) - lfus )
+                                   ! but then qemp_ice should also include sublimation 
+      END DO
+
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
 #endif

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -91 +91 @@
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
+   REAL(wp), PUBLIC :: rn_sfac ! multiplication factor for snow precipitation over sea-ice
 
    !! * Substitutions
@@ -151 +152 @@
          &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
          &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
-         &                  sn_tdif, rn_zqt,  rn_zu
+         &                  sn_tdif, rn_zqt,  rn_zu, rn_sfac
       !!---------------------------------------------------------------------
       !
@@ -158 +159 @@
          !                                      ! ====================== !
          !
+         rn_sfac = 1._wp       ! Default to one if missing from namelist 
          REWIND( numnam_ref )              ! Namelist namsbc_core in reference namelist : CORE bulk parameters
          READ  ( numnam_ref, namsbc_core, IOSTAT = ios, ERR = 901)
@@ -206 +208 @@
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
          qlw_ice(:,:,1)   = sf(jp_qlw)%fnow(:,:,1) 
-         qsr_ice(:,:,1)   = sf(jp_qsr)%fnow(:,:,1)
+         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
+         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
+         ENDIF
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)         
          qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
@@ -403 +407 @@
          CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
          CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
+         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac   ! output total precipitation [kg/m2/s]
+         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac   ! output solid precipitation [kg/m2/s]
+         CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow
+         CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation
       ENDIF
       !
@@ -608 +616 @@
       ! --- evaporation --- !
       z1_lsub = 1._wp / Lsub
-      evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation
-      devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub
-      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean
+      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation
+      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT
+      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean
 
       ! --- evaporation minus precipitation --- !
@@ -633 +641 @@
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+
+      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
+      DO jl = 1, jpl
+         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) )
+                                   ! But we do not have Tice => consider it at 0°C => evap=0 
+      END DO
 
       CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -33 +33 @@
    USE cpl_oasis3      ! OASIS3 coupling
    USE geo2ocean       ! 
-   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
+   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev,            &
+                      CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl,            & 
+                      PCO2a_in_cpl, Dust_in_cpl, &
+                      ln_medusa
    USE albedo          !
    USE in_out_manager  ! I/O manager
@@ -46 +49 @@
    USE p4zflx, ONLY : oce_co2
 #endif
-#if defined key_cice
-   USE ice_domain_size, only: ncat
-#endif
 #if defined key_lim3
    USE limthd_dh       ! for CALL lim_thd_snwblow
 #endif
+   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -105 +106 @@
    INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness 
    INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_ts_ice = 43            ! skin temperature of sea-ice (used for melt-ponds)
+   INTEGER, PARAMETER ::   jpr_grnm   = 44            ! Greenland ice mass
+   INTEGER, PARAMETER ::   jpr_antm   = 45            ! Antarctic ice mass
+   INTEGER, PARAMETER ::   jpr_atm_pco2 = 46          ! Incoming atm CO2 flux
+   INTEGER, PARAMETER ::   jpr_atm_dust = 47          ! Incoming atm aggregate dust 
+   INTEGER, PARAMETER ::   jprcv      = 47            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -135 +141 @@
    INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
    INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended
+   INTEGER, PARAMETER ::   jps_a_p    = 29            ! meltpond fraction  
+   INTEGER, PARAMETER ::   jps_ht_p   = 30            ! meltpond depth (m) 
+   INTEGER, PARAMETER ::   jps_kice   = 31            ! ice surface layer thermal conductivity
+   INTEGER, PARAMETER ::   jps_sstfrz = 32            ! sea-surface freezing temperature
+   INTEGER, PARAMETER ::   jps_fice1  = 33            ! first-order ice concentration (for time-travelling ice coupling)
+   INTEGER, PARAMETER ::   jps_bio_co2 = 34           ! MEDUSA air-sea CO2 flux
+   INTEGER, PARAMETER ::   jps_bio_dms = 35           ! MEDUSA DMS surface concentration
+   INTEGER, PARAMETER ::   jps_bio_chloro = 36        ! MEDUSA chlorophyll surface concentration
+   INTEGER, PARAMETER ::   jpsnd      = 36            ! total number of fields sent
+
+   REAL(wp), PARAMETER :: dms_unit_conv = 1.0e+6      ! Coversion factor to get outgong DMS in standard units for coupling
+                                                 ! i.e. specifically nmol/L (= umol/m3)
 
    !                                                         !!** namelist namsbc_cpl **
@@ -146 +163 @@
    END TYPE FLD_C
    ! Send to the atmosphere                           !
-   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                        
+   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2, sn_snd_cond, sn_snd_mpnd, sn_snd_sstfrz, sn_snd_thick1
+   TYPE(FLD_C) ::   sn_snd_bio_co2, sn_snd_bio_dms, sn_snd_bio_chloro                   
+
    ! Received from the atmosphere                     !
    TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                        
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm
+   TYPE(FLD_C) ::   sn_rcv_atm_pco2, sn_rcv_atm_dust                         
+
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -188 +209 @@
       ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
 #endif
-      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
+      !ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
+      ! Hardwire only two models as nn_cplmodel has not been read in
+      ! from the namelist yet.
+      ALLOCATE( xcplmask(jpi,jpj,0:2) , STAT=ierr(3) )   
       !
       sbc_cpl_alloc = MAXVAL( ierr )
@@ -216 +240 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      &
-         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      &
-         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   &
-         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask
+      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick , sn_snd_crt   , sn_snd_co2,     &
+         &                  sn_snd_cond, sn_snd_mpnd  , sn_snd_sstfrz, sn_snd_thick1,                 &
+         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     &
+         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
+         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
+         &                  ln_usecplmask, 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, rn_iceshelf_fluxes_tolerance
       !!---------------------------------------------------------------------
+
+      ! Add MEDUSA related fields to namelist
+      NAMELIST/namsbc_cpl/  sn_snd_bio_co2, sn_snd_bio_dms, sn_snd_bio_chloro,                        &
+         &                  sn_rcv_atm_pco2, sn_rcv_atm_dust
+
+      !!---------------------------------------------------------------------
+
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
@@ -245 +280 @@
       ENDIF
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
-         WRITE(numout,*)'  received fields (mutiple ice categogies)'
+         WRITE(numout,*)'  received fields (mutiple ice categories)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
          WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
@@ -258 +293 @@
          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,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
+         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), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
@@ -268 +307 @@
          WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
+         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), ')'
          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 & depth      = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat   ), ')'
+         WRITE(numout,*)'      sea surface freezing temp       = ', TRIM(sn_snd_sstfrz%cldes   ), ' (', TRIM(sn_snd_sstfrz%clcat   ), ')'
+
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         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
 
       !                                   ! allocate sbccpl arrays
-      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
+      !IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
      
       ! ================================ !
@@ -337 +390 @@
          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' ) 
@@ -383 +440 @@
       srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
       srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
-      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
+      srcv(jpr_ievp)%clname = 'OIceEvp'      ! evaporation over ice = sublimation
       srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation 
       srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
@@ -396 +453 @@
       CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
       END SELECT
-
+      !Set the number of categories for coupling of sublimation
+      IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl
+      !
       !                                                      ! ------------------------- !
       !                                                      !     Runoffs & Calving     !   
@@ -410 +469 @@
       !
       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' )   srcv(jpr_grnm)%laction = .TRUE.
+      srcv(jpr_antm  )%clname = 'OAntmass'   ;   IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' )   srcv(jpr_antm)%laction = .TRUE.
+
 
       !                                                      ! ------------------------- !
@@ -470 +532 @@
       !                                                      ! ------------------------- !
       srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
+
+
+      !                                                      ! --------------------------------------- !    
+      !                                                      ! 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
+    
       !                                                      ! ------------------------- !
       !                                                      !   topmelt and botmelt     !   
@@ -483 +561 @@
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
+      
+#if defined key_cice && ! defined key_cice4
+      !                                                      ! ----------------------------- !
+      !                                                      !  sea-ice skin temperature     !   
+      !                                                      !  used in meltpond scheme      !
+      !                                                      !  May be calculated in Atm     !
+      !                                                      ! ----------------------------- !
+      srcv(jpr_ts_ice)%clname = 'OTsfIce'
+      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' ) srcv(jpr_ts_ice)%laction = .TRUE.
+      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' ) srcv(jpr_ts_ice)%nct = jpl
+      !TODO: Should there be a consistency check here?
+#endif
+
       !                                                      ! ------------------------------- !
       !                                                      !   OPA-SAS coupling - rcv by opa !   
@@ -600 +691 @@
       !                                                      ! ------------------------- !
       ssnd(jps_toce)%clname = 'O_SSTSST'
-      ssnd(jps_tice)%clname = 'O_TepIce'
+      ssnd(jps_tice)%clname = 'OTepIce'
       ssnd(jps_tmix)%clname = 'O_TepMix'
       SELECT CASE( TRIM( sn_snd_temp%cldes ) )
       CASE( 'none'                                 )       ! nothing to do
       CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
-      CASE( 'oce and ice' , 'weighted oce and ice' )
+      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice')
          ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
          IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
@@ -634 +725 @@
 
       !                                                      ! ------------------------- !
-      !                                                      !  Ice fraction & Thickness ! 
+      !                                                      !  Ice fraction & Thickness 
       !                                                      ! ------------------------- !
       ssnd(jps_fice)%clname = 'OIceFrc'
       ssnd(jps_hice)%clname = 'OIceTck'
       ssnd(jps_hsnw)%clname = 'OSnwTck'
+      ssnd(jps_a_p)%clname  = 'OPndFrc'
+      ssnd(jps_ht_p)%clname = 'OPndTck'
+      ssnd(jps_fice1)%clname = 'OIceFrd'
       IF( k_ice /= 0 ) THEN
          ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
+         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used
+                                                     ! in producing atmos-to-ice fluxes
 ! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
          IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
+         IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = jpl
       ENDIF
       
@@ -657 +754 @@
       CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
       END SELECT
+
+      !                                                      ! ------------------------- !
+      !                                                      ! Ice Meltponds             !
+      !                                                      ! ------------------------- !
+#if defined key_cice && ! defined key_cice4
+      ! Meltponds only CICE5 
+      ssnd(jps_a_p)%clname = 'OPndFrc'   
+      ssnd(jps_ht_p)%clname = 'OPndTck'   
+      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) )
+      CASE ( 'none' )
+         ssnd(jps_a_p)%laction = .FALSE.
+         ssnd(jps_ht_p)%laction = .FALSE.
+      CASE ( 'ice only' ) 
+         ssnd(jps_a_p)%laction = .TRUE.
+         ssnd(jps_ht_p)%laction = .TRUE.
+         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
+            ssnd(jps_a_p)%nct = jpl
+            ssnd(jps_ht_p)%nct = jpl
+         ELSE
+            IF ( jpl > 1 ) THEN
+               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' )
+            ENDIF
+         ENDIF
+      CASE ( 'weighted ice' ) 
+         ssnd(jps_a_p)%laction = .TRUE.
+         ssnd(jps_ht_p)%laction = .TRUE.
+         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
+            ssnd(jps_a_p)%nct = jpl 
+            ssnd(jps_ht_p)%nct = jpl 
+         ENDIF
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes' )
+      END SELECT
+#else
+      IF( TRIM( sn_snd_mpnd%cldes ) /= 'none' ) THEN
+         CALL ctl_stop('Meltponds can only be used with CICEv5')
+      ENDIF
+#endif
 
       !                                                      ! ------------------------- !
@@ -689 +823 @@
       !                                                      ! ------------------------- !
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
+      !
+
+      !                                                      ! ------------------------- !
+      !                                                      !   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.
+
+      !                                                      ! ------------------------- !
+      !                                                      ! Sea surface freezing temp !
+      !                                                      ! ------------------------- !
+      ssnd(jps_sstfrz)%clname = 'O_SSTFrz' ; IF( TRIM(sn_snd_sstfrz%cldes) == 'coupled' )  ssnd(jps_sstfrz)%laction = .TRUE.
+      !
+      !                                                      ! ------------------------- !
+      !                                                      !    Ice conductivity       !
+      !                                                      ! ------------------------- !
+      ! Note that ultimately we will move to passing an ocean effective conductivity as well so there
+      ! will be some changes to the parts of the code which currently relate only to ice conductivity
+      ssnd(jps_kice )%clname = 'OIceKn'
+      SELECT CASE ( TRIM( sn_snd_cond%cldes ) )
+      CASE ( 'none' )
+         ssnd(jps_kice)%laction = .FALSE.
+      CASE ( 'ice only' )
+         ssnd(jps_kice)%laction = .TRUE.
+         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN
+            ssnd(jps_kice)%nct = jpl
+         ELSE
+            IF ( jpl > 1 ) THEN
+               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' )
+            ENDIF
+         ENDIF
+      CASE ( 'weighted ice' )
+         ssnd(jps_kice)%laction = .TRUE.
+         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = jpl
+      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes' )
+      END SELECT
+      !
+      
 
       !                                                      ! ------------------------------- !
@@ -785 +966 @@
       ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
+      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
+
+          ! initialise other variables
+          greenland_icesheet_mass_array(:,:) = 0.0
+          antarctica_icesheet_mass_array(:,:) = 0.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
+
       CALL wrk_dealloc( jpi,jpj, zacs, zaos )
       !
@@ -843 +1047 @@
       !!
       LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
-      INTEGER  ::   ji, jj, jn             ! dummy loop indices
+      INTEGER  ::   ji, jj, jl, jn         ! dummy loop indices
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
+      INTEGER  ::   ikchoix
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
+      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) ::   zcoef                  ! temporary scalar
       REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
       REAL(wp) ::   zzx, zzy               ! temporary variables
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2
       !!----------------------------------------------------------------------
+
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
       !
-      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
+      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2 )
       !
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -893 +1102 @@
             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 )   
-               ELSE  
-                  CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%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
+             DO jj=2,jpjm1
+                DO ji=2,jpim1
+                     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))
+                ENDDO
+             ENDDO
+                   
+             ikchoix = 1
+             CALL repcmo (frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)
+             CALL lbc_lnk (ztx2,'U', -1. )
+             CALL lbc_lnk (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->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
             !                              
@@ -990 +1221 @@
       ENDIF
 
+      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
+
 #if defined key_cpl_carbon_cycle
       !                                                      ! ================== !
@@ -995 +1231 @@
       !                                                      ! ================== !
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
+#endif
+
+#if defined key_cice && ! defined key_cice4
+      !  ! Sea ice surface skin temp:
+      IF( srcv(jpr_ts_ice)%laction ) THEN
+        DO jl = 1, jpl
+          DO jj = 1, jpj
+            DO ji = 1, jpi
+              IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) > 0.0) THEN
+                tsfc_ice(ji,jj,jl) = 0.0
+              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) < -60.0) THEN
+                tsfc_ice(ji,jj,jl) = -60.0
+              ELSE
+                tsfc_ice(ji,jj,jl) = frcv(jpr_ts_ice)%z3(ji,jj,jl)
+              ENDIF
+            END DO
+          END DO
+        END DO
+      ENDIF
 #endif
 
@@ -1029 +1284 @@
          ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
          ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
+         un (:,:,1) = ssu_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
          CALL iom_put( 'ssu_m', ssu_m )
       ENDIF
@@ -1034 +1290 @@
          ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
          vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
+         vn (:,:,1) = ssv_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
          CALL iom_put( 'ssv_m', ssv_m )
       ENDIF
@@ -1110 +1367 @@
 
       ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
+      
+      !                                                        ! land ice masses : Greenland
+      zepsilon = rn_iceshelf_fluxes_tolerance
+
+
+      ! See if we need zmask_sum...
+      IF ( srcv(jpr_grnm)%laction .OR. srcv(jpr_antm)%laction ) THEN
+         zmask_sum = glob_sum( tmask(:,:,1) )
+      ENDIF
+
+      IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
+         greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
+         ! take average over ocean points of input array to avoid cumulative error over time
+         ! The following must be bit reproducible over different PE decompositions
+         zgreenland_icesheet_mass_in = glob_sum( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
+
+         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum
+         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) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
+      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
+         antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
+         ! take average over ocean points of input array to avoid cumulative error from rounding errors over time
+         ! The following must be bit reproducible over different PE decompositions
+         zantarctica_icesheet_mass_in = glob_sum( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
+
+         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum
+         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) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
+      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
+         antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
+      ENDIF
+
+      !
+      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
@@ -1333 +1673 @@
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
       !!
-      !! ** Purpose :   provide the heat and freshwater fluxes of the 
-      !!              ocean-ice system.
+      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
       !!
       !! ** Method  :   transform the fields received from the atmosphere into
       !!             surface heat and fresh water boundary condition for the 
       !!             ice-ocean system. The following fields are provided:
-      !!              * total non solar, solar and freshwater fluxes (qns_tot, 
+      !!               * total non solar, solar and freshwater fluxes (qns_tot, 
       !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
       !!             NB: emp_tot include runoffs and calving.
-      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
+      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
       !!             emp_ice = sublimation - solid precipitation as liquid
       !!             precipitation are re-routed directly to the ocean and 
-      !!             runoffs and calving directly enter the ocean.
-      !!              * solid precipitation (sprecip), used to add to qns_tot 
+      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
+      !!               * solid precipitation (sprecip), used to add to qns_tot 
       !!             the heat lost associated to melting solid precipitation
       !!             over the ocean fraction.
-      !!       ===>> CAUTION here this changes the net heat flux received from
-      !!             the atmosphere
-      !!
-      !!                  - the fluxes have been separated from the stress as
-      !!                 (a) they are updated at each ice time step compare to
-      !!                 an update at each coupled time step for the stress, and
-      !!                 (b) the conservative computation of the fluxes over the
-      !!                 sea-ice area requires the knowledge of the ice fraction
-      !!                 after the ice advection and before the ice thermodynamics,
-      !!                 so that the stress is updated before the ice dynamics
-      !!                 while the fluxes are updated after it.
+      !!               * heat content of rain, snow and evap can also be provided,
+      !!             otherwise heat flux associated with these mass flux are
+      !!             guessed (qemp_oce, qemp_ice)
+      !!
+      !!             - the fluxes have been separated from the stress as
+      !!               (a) they are updated at each ice time step compare to
+      !!               an update at each coupled time step for the stress, and
+      !!               (b) the conservative computation of the fluxes over the
+      !!               sea-ice area requires the knowledge of the ice fraction
+      !!               after the ice advection and before the ice thermodynamics,
+      !!               so that the stress is updated before the ice dynamics
+      !!               while the fluxes are updated after it.
+      !!
+      !! ** Details
+      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
+      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
+      !!
+      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
+      !!
+      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce)
+      !!                                                                      river runoff (rnf) is provided but not included here
       !!
       !! ** Action  :   update at each nf_ice time step:
       !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
       !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
-      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
-      !!                   emp_ice            ice sublimation - solid precipitation over the ice
-      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
-      !!                   sprecip             solid precipitation over the ocean  
+      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
+      !!                   emp_ice           ice sublimation - solid precipitation over the ice
+      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
+      !!                   sprecip           solid precipitation over the ocean  
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
@@ -1376 +1725 @@
       !
       INTEGER ::   jl         ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk, zsnw
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
+      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
+      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
 
       IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
@@ -1392 +1743 @@
       !
       !                                                      ! ========================= !
-      !                                                      !    freshwater budget      !   (emp)
+      !                                                      !    freshwater budget      !   (emp_tot)
       !                                                      ! ========================= !
       !
-      !                                                           ! total Precipitation - total Evaporation (emp_tot)
-      !                                                           ! solid precipitation - sublimation       (emp_ice)
-      !                                                           ! solid Precipitation                     (sprecip)
-      !                                                           ! liquid + solid Precipitation            (tprecip)
+      !                                                           ! solid Precipitation                                (sprecip)
+      !                                                           ! liquid + solid Precipitation                       (tprecip)
+      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
+      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
       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(:,:)
-         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
-            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation 
+         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)         
+#if defined key_cice
+         IF ( TRIM(sn_rcv_emp%clcat) == 'yes' ) THEN
+            ! zemp_ice is the sum of frcv(jpr_ievp)%z3(:,:,1) over all layers - snow
+            zemp_ice(:,:) = - frcv(jpr_snow)%z3(:,:,1)
+            DO jl=1,jpl
+               zemp_ice(:,:   ) = zemp_ice(:,:) + frcv(jpr_ievp)%z3(:,:,jl)
+            ENDDO
+            ! latent heat coupled for each category in CICE
+            qla_ice(:,:,1:jpl) = - frcv(jpr_ievp)%z3(:,:,1:jpl) * lsub
+         ELSE
+            ! If CICE has multicategories it still expects coupling fields for
+            ! each even if we treat as a single field
+            ! The latent heat flux is split between the ice categories according
+            ! to the fraction of the ice in each category
+            zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
+            WHERE ( zicefr(:,:) /= 0._wp ) 
+               ztmp(:,:) = 1./zicefr(:,:)
+            ELSEWHERE 
+               ztmp(:,:) = 0.e0
+            END WHERE  
+            DO jl=1,jpl
+               qla_ice(:,:,jl) = - a_i(:,:,jl) * ztmp(:,:) * frcv(jpr_ievp)%z3(:,:,1) * lsub 
+            END DO
+            WHERE ( zicefr(:,:) == 0._wp )  qla_ice(:,:,1) = -frcv(jpr_ievp)%z3(:,:,1) * lsub 
+         ENDIF
+#else         
+         zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * zicefr(:,:)
+#endif                  
+         CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1) * tmask(:,:,1)      )   ! liquid precipitation 
+         CALL iom_put( 'rain_ao_cea'  , frcv(jpr_rain)%z3(:,:,1)* p_frld(:,:) * tmask(:,:,1)      )   ! liquid precipitation 
          IF( iom_use('hflx_rain_cea') )   &
-            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip. 
+            &  CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) * tmask(:,:,1))   ! heat flux from liq. precip. 
+         IF( iom_use('hflx_prec_cea') )   &
+            & CALL iom_put( 'hflx_prec_cea', ztprecip * zcptn(:,:) * tmask(:,:,1) * p_frld(:,:) )   ! heat content flux from all precip  (cell avg)
          IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
-            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+            & ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
          IF( iom_use('evap_ao_cea'  ) )   &
-            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
+            &  CALL iom_put( 'evap_ao_cea'  , ztmp * tmask(:,:,1)                  )   ! ice-free oce evap (cell average)
          IF( iom_use('hflx_evap_cea') )   &
-            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
-      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
+            &  CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) * tmask(:,:,1) )   ! heat flux from from evap (cell average)
+      CASE( 'oce and ice' )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
          zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
-         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
+         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * zicefr(:,:)
          zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
          ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
       END SELECT
 
-      IF( iom_use('subl_ai_cea') )   &
-         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
-      !   
-      !                                                           ! runoffs and calving (put in emp_tot)
+#if defined key_lim3
+      ! zsnw = snow fraction over ice after wind blowing
+      zsnw(:,:) = 0._wp  ;  CALL lim_thd_snwblow( p_frld, zsnw )
+      
+      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
+      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( zicefr(:,:) - 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) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
+
+      ! --- evaporation over ice (kg/m2/s) --- !
+      zevap_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1)
+      ! 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 in case no subgrid scale fluxes are provided by atm.
+      zdevap_ice(:,:) = 0._wp
+      
+      ! --- runoffs (included in emp later on) --- !
       IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+
+      ! --- calving (put in emp_tot and emp_oce) --- !
+      IF( srcv(jpr_cal)%laction ) THEN 
+         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
+      ENDIF
+
+      IF( ln_mixcpl ) THEN
+         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
+         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
+         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:)
+         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
+         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
+         DO jl=1,jpl
+            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:)
+            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:)
+         ENDDO
+      ELSE
+         emp_tot(:,:) =         zemp_tot(:,:)
+         emp_ice(:,:) =         zemp_ice(:,:)
+         emp_oce(:,:) =         zemp_oce(:,:)     
+         sprecip(:,:) =         zsprecip(:,:)
+         tprecip(:,:) =         ztprecip(:,:)
+         DO jl=1,jpl
+            evap_ice (:,:,jl) = zevap_ice (:,:)
+            devap_ice(:,:,jl) = zdevap_ice(:,:)
+         ENDDO
+      ENDIF
+
+      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea', zevap_ice(:,:) * zicefr(:,:)         )  ! Sublimation over sea-ice (cell average)
+                                     CALL iom_put( 'snowpre'    , sprecip(:,:)                         )  ! Snow
+      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea', sprecip(:,:) * ( 1._wp - zsnw(:,:) ) )  ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea', sprecip(:,:) *           zsnw(:,:)   )  ! Snow over sea-ice         (cell average)
+#else
+      ! runoffs and calving (put in emp_tot)
+      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      IF( iom_use('hflx_rnf_cea') )   &
+         CALL iom_put( 'hflx_rnf_cea' , rnf(:,:) * zcptn(:,:) )
       IF( srcv(jpr_cal)%laction ) THEN 
          zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
@@ -1443 +1877 @@
       ENDIF
 
-         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
-      IF( iom_use('snow_ao_cea') )   &
-         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
-      IF( iom_use('snow_ai_cea') )   &
-         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
+      IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )  ! Sublimation over sea-ice (cell average)
+                                    CALL iom_put( 'snowpre'    , sprecip(:,:)               )   ! Snow
+      IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:) )   ! Snow over ice-free ocean  (cell average)
+      IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:) )   ! Snow over sea-ice         (cell average)
+#endif
 
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
       !                                                      ! ========================= !
-      CASE( 'oce only' )                                     ! the required field is directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
-      CASE( 'conservative' )                                      ! the required fields are directly provided
-         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
+      CASE( 'oce only' )         ! the required field is directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'conservative' )     ! the required fields are directly provided
+         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
-            ! Set all category values equal for the moment
             DO jl=1,jpl
-               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
             ENDDO
          ENDIF
-      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
-         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
+      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
+         zqns_tot(:,:) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
@@ -1472 +1905 @@
             ENDDO
          ELSE
-            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            qns_tot(:,:) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
                zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
@@ -1478 +1911 @@
             ENDDO
          ENDIF
-      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
+      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
          zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
          zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
             &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
-            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
+            &                                           + pist(:,:,1) * zicefr(:,:) ) )
       END SELECT
 !!gm
@@ -1493 +1926 @@
 !! similar job should be done for snow and precipitation temperature
       !                                     
-      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting 
-         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting 
-         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
-         IF( iom_use('hflx_cal_cea') )   &
-            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
-      ENDIF
-
-      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
-      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
-
-#if defined key_lim3
-      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
-
-      ! --- evaporation --- !
-      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation
-      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice
-      !                 but it is incoherent WITH the ice model  
-      DO jl=1,jpl
-         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1)
-      ENDDO
-      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
-
-      ! --- evaporation minus precipitation --- !
-      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:)
-
+      IF( srcv(jpr_cal)%laction ) THEN   ! Iceberg melting 
+         zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! add the latent heat of iceberg melting
+                                                                         ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
+         IF( iom_use('hflx_cal_cea') )   CALL iom_put( 'hflx_cal_cea', - frcv(jpr_cal)%z3(:,:,1) * lfus )   ! heat flux from calving
+      ENDIF
+
+#if defined key_lim3      
       ! --- non solar flux over ocean --- !
       !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
@@ -1523 +1938 @@
       WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
 
-      ! --- heat flux associated with emp --- !
-      zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
-      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
-         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
-         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
-      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
-         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
-
-      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
+      ! --- heat flux associated with emp (W/m2) --- !
+      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn(:,:)   &       ! evap
+         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &       ! liquid precip
+         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus )  ! solid precip over ocean + snow melting
+!      zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
+!         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
+      zqemp_ice(:,:) =      zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice (only)
+                                                                                                       ! qevap_ice=0 since we consider Tice=0degC
+      
+      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
       zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
 
-      ! --- total non solar flux --- !
-      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
+      ! --- 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 ) * cpic ) but we do not have Tice, so we consider Tice=0degC
+      END DO
+
+      ! --- total non solar flux (including evap/precip) --- !
+      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)
 
       ! --- in case both coupled/forced are active, we must mix values --- ! 
@@ -1543 +1963 @@
          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(:,:)
+            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
+            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
          ENDDO
          qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
          qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
-!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0)
+         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
       ELSE
          qns_tot  (:,:  ) = zqns_tot  (:,:  )
          qns_oce  (:,:  ) = zqns_oce  (:,:  )
          qns_ice  (:,:,:) = zqns_ice  (:,:,:)
-         qprec_ice(:,:)   = zqprec_ice(:,:)
-         qemp_oce (:,:)   = zqemp_oce (:,:)
-      ENDIF
-
-      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
+         qevap_ice(:,:,:) = zqevap_ice(:,:,:)
+         qprec_ice(:,:  ) = zqprec_ice(:,:  )
+         qemp_oce (:,:  ) = zqemp_oce (:,:  )
+         qemp_ice (:,:  ) = zqemp_ice (:,:  )
+      ENDIF
+
+      !! clem: we should output qemp_oce and qemp_ice (at least)
+      IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average)
+      !! these diags are not outputed yet
+!!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average)
+!!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average)
+!!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average)
+
 #else
-
       ! clem: this formulation is certainly wrong... but better than it was...
+      
       zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
-         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
+         &          - (p_frld(:,:) * zsprecip(:,:) * lfus)  &          ! remove the latent heat flux of solid precip. melting
          &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
-         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
+         &             - zemp_ice(:,:) ) * zcptn(:,:) 
 
      IF( ln_mixcpl ) THEN
@@ -1575 +2004 @@
          qns_ice(:,:,:) = zqns_ice(:,:,:)
       ENDIF
-
 #endif
 
@@ -1626 +2054 @@
 
 #if defined key_lim3
-      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
       ! --- solar flux over ocean --- !
       !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
@@ -1634 +2061 @@
       IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
       ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
-
-      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
 #endif
 
@@ -1686 +2111 @@
       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 
-      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
-      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
+      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zsnw )
+      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice )
+      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice )
+      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
@@ -1706 +2133 @@
       !
       INTEGER ::   ji, jj, jl   ! dummy loop indices
+      INTEGER ::   ikchoix
       INTEGER ::   isec, info   ! local integer
       REAL(wp) ::   zumax, zvmax
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zotx1_in, zoty1_in
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
       !!----------------------------------------------------------------------
@@ -1715 +2144 @@
       !
       CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
+      CALL wrk_alloc( jpi,jpj, zotx1_in, zoty1_in)
       CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
 
@@ -1743 +2173 @@
                      ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
                   ELSEWHERE
-                     ztmp3(:,:,1) = rt0 ! TODO: Is freezing point a good default? (Maybe SST is better?)
+                     ztmp3(:,:,1) = rt0
                   END WHERE
                CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
@@ -1758 +2188 @@
                CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
                END SELECT
+            CASE( 'oce and weighted ice' )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0 
+               SELECT CASE( sn_snd_temp%clcat )
+               CASE( 'yes' )   
+           ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
+               CASE( 'no' )
+           ztmp3(:,:,:) = 0.0
+           DO jl=1,jpl
+                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
+           ENDDO
+               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
+               END SELECT
             CASE( 'mixed oce-ice'        )   
                ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
@@ -1774 +2215 @@
       !                                                      ! ------------------------- !
       IF( ssnd(jps_albice)%laction ) THEN                         ! ice 
-         SELECT CASE( sn_snd_alb%cldes )
-         CASE( 'ice'          )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
-         CASE( 'weighted ice' )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
-         CASE default             ; CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
+          SELECT CASE( sn_snd_alb%cldes )
+          CASE( 'ice' )
+             SELECT CASE( sn_snd_alb%clcat )
+             CASE( 'yes' )   
+                ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
+             CASE( 'no' )
+                WHERE( SUM( a_i, dim=3 ) /= 0. )
+                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
+                ELSEWHERE
+                   ztmp1(:,:) = albedo_oce_mix(:,:)
+                END WHERE
+             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
+             END SELECT
+          CASE( 'weighted ice' )   ;
+             SELECT CASE( sn_snd_alb%clcat )
+             CASE( 'yes' )   
+                ztmp3(:,:,1:jpl) =  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
+             CASE( 'no' )
+                WHERE( fr_i (:,:) > 0. )
+                   ztmp1(:,:) = SUM (  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
+                ELSEWHERE
+                   ztmp1(:,:) = 0.
+                END WHERE
+             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ice%clcat' )
+             END SELECT
+          CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
          END SELECT
-         CALL cpl_snd( jps_albice, isec, ztmp3, info )
-      ENDIF
+
+         SELECT CASE( sn_snd_alb%clcat )
+            CASE( 'yes' )   
+               CALL cpl_snd( jps_albice, isec, ztmp3, info )      !-> MV this has never been checked in coupled mode
+            CASE( 'no'  )   
+               CALL cpl_snd( jps_albice, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
+         END SELECT
+      ENDIF
+
       IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
          ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
@@ -1799 +2269 @@
          END SELECT
          IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
+      ENDIF
+      
+      ! Send ice fraction field (first order interpolation), for weighting UM fluxes to be passed to NEMO
+      IF (ssnd(jps_fice1)%laction) THEN
+         SELECT CASE (sn_snd_thick1%clcat)
+         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
+         CASE( 'no' )    ;   ztmp3(:,:,1) = fr_i(:,:)
+         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' )
+    END SELECT
+         CALL cpl_snd (jps_fice1, isec, ztmp3, info)
       ENDIF
       
@@ -1845 +2325 @@
       ENDIF
       !
+#if defined key_cice && ! defined key_cice4
+      ! Send meltpond fields 
+      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN
+         SELECT CASE( sn_snd_mpnd%cldes) 
+         CASE( 'weighted ice' ) 
+            SELECT CASE( sn_snd_mpnd%clcat ) 
+            CASE( 'yes' ) 
+               ztmp3(:,:,1:jpl) =  a_p(:,:,1:jpl) * a_i(:,:,1:jpl) 
+               ztmp4(:,:,1:jpl) =  ht_p(:,:,1:jpl) * a_i(:,:,1:jpl) 
+            CASE( 'no' ) 
+               ztmp3(:,:,:) = 0.0 
+               ztmp4(:,:,:) = 0.0 
+               DO jl=1,jpl 
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_p(:,:,jpl) * a_i(:,:,jpl) 
+                 ztmp4(:,:,1) = ztmp4(:,:,1) + ht_p(:,:,jpl) * a_i(:,:,jpl) 
+               ENDDO 
+            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' ) 
+            END SELECT 
+         CASE( 'ice only' )    
+            ztmp3(:,:,1:jpl) = a_p(:,:,1:jpl) 
+            ztmp4(:,:,1:jpl) = ht_p(:,:,1:jpl) 
+         END SELECT 
+         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )    
+         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )    
+         !
+         ! Send ice effective conductivity
+         SELECT CASE( sn_snd_cond%cldes)
+         CASE( 'weighted ice' )   
+            SELECT CASE( sn_snd_cond%clcat )
+            CASE( 'yes' )   
+               ztmp3(:,:,1:jpl) =  kn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
+            CASE( 'no' )
+               ztmp3(:,:,:) = 0.0
+               DO jl=1,jpl
+                 ztmp3(:,:,1) = ztmp3(:,:,1) + kn_ice(:,:,jl) * a_i(:,:,jl)
+               ENDDO
+            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' )
+            END SELECT
+         CASE( 'ice only' )   
+           ztmp3(:,:,1:jpl) = kn_ice(:,:,1:jpl)
+         END SELECT
+         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info )
+      ENDIF
+#endif
+      !
+      !
 #if defined key_cpl_carbon_cycle
       !                                                      ! ------------------------- !
@@ -1852 +2378 @@
       !
 #endif
+
+
+
+      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
+
       !                                                      ! ------------------------- !
       IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
@@ -1858 +2404 @@
          !                                                  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-------|
@@ -1867 +2413 @@
             zotx1(:,:) = un(:,:,1)  
             zoty1(:,:) = vn(:,:,1)  
-         ELSE        
+         ELSE
             SELECT CASE( TRIM( sn_snd_crt%cldes ) )
             CASE( 'oce only'             )      ! C-grid ==> T
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
-                     zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
+               IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
+                  DO jj = 2, jpjm1
+                     DO ji = fs_2, fs_jpim1   ! vector opt.
+                        zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
+                        zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
+                     END DO
                   END DO
-               END DO
+               ELSE
+! Temporarily Changed for UKV
+                  DO jj = 2, jpjm1
+                     DO ji = 2, jpim1
+                        zotx1(ji,jj) = un(ji,jj,1)
+                        zoty1(ji,jj) = vn(ji,jj,1)
+                     END DO
+                  END DO
+               ENDIF 
             CASE( 'weighted oce and ice' )   
                SELECT CASE ( cp_ice_msh )
@@ -1934 +2490 @@
                   END DO
                CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
-                  DO jj = 2, jpjm1
-                     DO ji = 2, jpim1   ! NO vector opt.
-                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
-                           &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
-                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
-                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
-                           &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
-                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
+                  IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
+                     DO jj = 2, jpjm1
+                        DO ji = 2, jpim1   ! NO vector opt.
+                           zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj,1) ) * zfr_l(ji,jj)   &   
+                                &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
+                                &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
+                           zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji,jj-1,1) ) * zfr_l(ji,jj)   &
+                                &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
+                                &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
+                        END DO
                      END DO
-                  END DO
+#if defined key_cice
+                  ELSE
+! Temporarily Changed for HadGEM3
+                     DO jj = 2, jpjm1
+                        DO ji = 2, jpim1   ! NO vector opt.
+                           zotx1(ji,jj) = (1.0-fr_iu(ji,jj)) * un(ji,jj,1)             &
+                                &              + fr_iu(ji,jj) * 0.5 * ( u_ice(ji,jj-1) + u_ice(ji,jj) ) 
+                           zoty1(ji,jj) = (1.0-fr_iv(ji,jj)) * vn(ji,jj,1)             &
+                                &              + fr_iv(ji,jj) * 0.5 * ( v_ice(ji-1,jj) + v_ice(ji,jj) ) 
+                        END DO
+                     END DO
+#endif
+                  ENDIF
                END SELECT
             END SELECT
@@ -1953 +2523 @@
          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(:,:)
-            ENDIF
+            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 jj=2,jpjm1
+                  DO ji=2,jpim1
+                     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))
+                  ENDDO
+               ENDDO
+
+               ! Ensure any N fold and wrap columns are updated
+               CALL lbc_lnk(ztmp1, 'V', -1.0)
+               CALL lbc_lnk(ztmp2, 'U', -1.0)
+                            
+               ikchoix = -1
+               ! We need copies of zotx1 and zoty2 in order to avoid problems 
+               ! caused by INTENTs used in the following subroutine. 
+               zotx1_in(:,:) = zotx1(:,:)
+               zoty1_in(:,:) = zoty1(:,:)
+               CALL repcmo (zotx1_in,ztmp2,ztmp1,zoty1_in,zotx1,zoty1,ikchoix)
+           ENDIF
          ENDIF
          !
@@ -2023 +2621 @@
       IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
       IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
-
+      
+#if defined key_cice
+      ztmp1(:,:) = sstfrz(:,:) + rt0
+      IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
+#endif
+      !
       CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
+      CALL wrk_dealloc( jpi,jpj, zotx1_in, zoty1_in )
       CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
       !

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

@@ -108 +108 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area   ! sum over the global domain
             zcoef = z_fwf * rcp
             emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1)
@@ -162 +162 @@
             zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
             !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges) 
-            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
+            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) - snwice_fmass(:,:) ) ) / area
             !            
             IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

@@ -15 +15 @@
    USE dom_oce         ! ocean space and time domain
    USE domvvl
-   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
+   USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
+   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
    USE in_out_manager  ! I/O manager
    USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit
@@ -37 +38 @@
    USE ice_gather_scatter
    USE ice_calendar, only: dt
+# if defined key_cice4
    USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
-# if defined key_cice4
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
                 strocnxT,strocnyT,                               & 
@@ -45 +46 @@
                 flatn_f,fsurfn_f,fcondtopn_f,                    &
                 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
-                swvdr,swvdf,swidr,swidf
+                swvdr,swvdf,swidr,swidf,Tf
    USE ice_therm_vertical, only: calc_Tsfc
 #else
+   USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
+                vsnon,vice,vicen,nt_Tsfc
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
                 strocnxT,strocnyT,                               & 
-                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,     &
-                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,          &
+                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,      &
+                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,             &
                 flatn_f,fsurfn_f,fcondtopn_f,                    &
+#ifdef key_asminc
+                daice_da,fresh_da,fsalt_da,                    &
+#endif
                 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
-                swvdr,swvdf,swidr,swidf
-   USE ice_therm_shared, only: calc_Tsfc
+                swvdr,swvdf,swidr,swidf,Tf,                      &
+      !! When using NEMO with CICE, this change requires use of 
+      !! one of the following two CICE branches:
+      !! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7
+      !! - at CICE5.1.2, hadax/vn5.1.2_GSI8
+                keffn_top,Tn_top
+
+   USE ice_therm_shared, only: calc_Tsfc, heat_capacity
+   USE ice_shortwave, only: apeffn
 #endif
    USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
@@ -161 +174 @@
       INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d
       INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices
       !!---------------------------------------------------------------------
@@ -174 +187 @@
       jj_off = INT ( (jpjglo - ny_global) / 2 )
 
-#if defined key_nemocice_decomp
-      ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
-      ! there is no restart file.
-      ! Values from a CICE restart file would overwrite this
-      IF ( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
-      ENDIF  
-#endif
-
-! Initialize CICE
+      ! Initialize CICE
       CALL CICE_Initialize
 
-! Do some CICE consistency checks
+      ! Do some CICE consistency checks
       IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          IF ( calc_strair .OR. calc_Tsfc ) THEN
@@ -198 +202 @@
 
 
-! allocate sbc_ice and sbc_cice arrays
-      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
+      ! allocate sbc_ice and sbc_cice arrays
+      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' )
       IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
 
-! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
+      ! Ensure that no temperature points are below freezing if not a NEMO restart
       IF( .NOT. ln_rstart ) THEN
-         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
+
+         CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d ) 
+         DO jk=1,jpk
+             CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), fsdept_n(:,:,jk) )
+         ENDDO
+         tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
          tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
-      ENDIF
-
-      fr_iu(:,:)=0.0
-      fr_iv(:,:)=0.0
+         CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d ) 
+
+#if defined key_nemocice_decomp
+         ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
+         ! there is no restart file.
+         ! Values from a CICE restart file would overwrite this
+         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
+#endif
+
+      ENDIF  
+
+      ! calculate surface freezing temperature and send to CICE
+      CALL  eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1)) 
+      CALL nemo2cice(sstfrz,Tf, 'T', 1. )
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
@@ -220 +239 @@
 ! T point to U point
 ! T point to V point
+      fr_iu(:,:)=0.0
+      fr_iv(:,:)=0.0
       DO jj=1,jpjm1
          DO ji=1,jpim1
@@ -283 +304 @@
  
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
-      !
+
+#if defined key_asminc
+      ! Initialize fresh water and salt fluxes from data assim   
+      !  and data assimilation index to cice 
+      nfresh_da(:,:) = 0.0   
+      nfsalt_da(:,:) = 0.0   
+      ndaice_da(:,:) = 0.0         
+#endif
+      !
+      ! In coupled mode get extra fields from CICE for passing back to atmosphere
+ 
+      IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)
+      ! 
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
       !
@@ -343 +376 @@
          CALL nemo2cice(ztmp,stray,'F', -1. )
 
+
+! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
+! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby 
+! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
+! gridbox mean fluxes in the UM by future ice concentration obtained through  
+! OASIS.  This allows for a much more realistic apportionment of energy through
+! the ice - and conserves energy.
+! Therefore the fluxes are now divided by ice concentration in the coupled
+! formulation (jp_purecpl) as well as for jp_flx.  This NEMO branch should only
+! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
+! which point the GSI8 UM changes were committed.
+
 ! Surface downward latent heat flux (CI_5)
          IF (ksbc == jp_flx) THEN
@@ -348 +393 @@
                ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
             ENDDO
-         ELSE
-! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
-            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
-! End of temporary code
-            DO jj=1,jpj
-               DO ji=1,jpi
-                  IF (fr_i(ji,jj).eq.0.0) THEN
-                     DO jl=1,ncat
-                        ztmpn(ji,jj,jl)=0.0
-                     ENDDO
-                     ! This will then be conserved in CICE
-                     ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
-                  ELSE
-                     DO jl=1,ncat
-                        ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
-                     ENDDO
-                  ENDIF
-               ENDDO
+         ELSE IF (ksbc == jp_purecpl) THEN
+            DO jl=1,ncat
+               ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
             ENDDO
+    ELSE
+           !In coupled mode - qla_ice calculated in sbc_cpl for each category
+           ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
          ENDIF
+
          DO jl=1,ncat
             CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
@@ -373 +407 @@
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
-            IF (ksbc == jp_flx) THEN
+            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
                ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
             ELSE
@@ -382 +416 @@
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
-            IF (ksbc == jp_flx) THEN
+            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
             ELSE
@@ -431 +465 @@
       ENDIF
 
+#if defined key_asminc
+!Ice concentration change (from assimilation)
+      ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
+      Call nemo2cice(ztmp,daice_da,'T', 1. )
+#endif 
+
 ! Snowfall
 ! 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)  
+      IF( kt == nit000 .AND. lwp )  THEN
+         WRITE(numout,*) 'sprecip weight, rn_sfac=', rn_sfac
+      ENDIF
+      ztmp(:,:)=MAX(fr_i(:,:)*rn_sfac*sprecip(:,:),0.0)  
       CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
 
@@ -442 +485 @@
       CALL nemo2cice(ztmp,frain,'T', 1. ) 
 
+! Recalculate freezing temperature and send to CICE 
+      CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1)) 
+      CALL nemo2cice(sstfrz,Tf,'T', 1. )
+
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
-      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
-
-      ztmp(:,:) = nfrzmlt(:,:)
-      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
+      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(sstfrz(:,:)-sst_m(:,:))/(2.0*dt) 
+      CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
 
 ! SST  and SSS
@@ -453 +498 @@
       CALL nemo2cice(sst_m,sst,'T', 1. )
       CALL nemo2cice(sss_m,sss,'T', 1. )
+
+      IF( ksbc == jp_purecpl ) THEN
+! Sea ice surface skin temperature
+         DO jl=1,ncat
+           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
+         ENDDO 
+      ENDIF
 
 ! x comp and y comp of surface ocean current
@@ -685 +737 @@
       ENDIF
 
+#if defined key_asminc
+! Import fresh water and salt flux due to seaice da
+      CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
+      CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
+#endif
+
 ! Release work space
 
@@ -708 +766 @@
       IF( nn_timing == 1 )  CALL timing_start('cice_sbc_hadgam')
       !
-      IF( kt == nit000 )  THEN
-         IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
-         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
-      ENDIF
-
       !                                         ! =========================== !
       !                                         !   Prepare Coupling fields   !
@@ -730 +783 @@
          CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
       ENDDO
+
+#if ! defined key_cice4
+! Meltpond fraction and depth
+      DO jl = 1,ncat
+         CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. )
+         CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. )
+      ENDDO
+#endif
+
+
+! If using multilayers thermodynamics in CICE then get top layer temperature
+! and effective conductivity       
+!! When using NEMO with CICE, this change requires use of 
+!! one of the following two CICE branches:
+!! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7
+!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
+      IF (heat_capacity) THEN
+         DO jl = 1,ncat
+            CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. )
+            CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. )
+         ENDDO
+! Convert surface temperature to Kelvin
+         tn_ice(:,:,:)=tn_ice(:,:,:)+rt0
+      ELSE
+         tn_ice(:,:,:) = 0.0
+         kn_ice(:,:,:) = 0.0
+      ENDIF       
+
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam')

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90

@@ -103 +103 @@
                                  ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
          
-         fr_i(:,:) = eos_fzp( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius]
+         CALL eos_fzp( sss_m(:,:), fr_i(:,:) )       ! sea surface freezing temperature [Celcius]
+         fr_i(:,:) = fr_i(:,:) * tmask(:,:,1)
 
          IF( ln_cpl )   a_i(:,:,1) = fr_i(:,:)         

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

@@ -110 +110 @@
       INTEGER  ::   jl                 ! dummy loop index
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled)
       REAL(wp), POINTER, DIMENSION(:,:  )   ::   zutau_ice, zvtau_ice 
       !!----------------------------------------------------------------------
@@ -126 +125 @@
          
          ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
-         t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )  
-         
+         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
+         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
+          
          ! Mask sea ice surface temperature (set to rt0 over land)
          DO jl = 1, jpl
@@ -196 +196 @@
          ! fr1_i0  , fr2_i0   : 1sr & 2nd fraction of qsr penetration in ice             [%]
          !----------------------------------------------------------------------------------------
-         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
          CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
 
@@ -202 +202 @@
          CASE( jp_clio )                                       ! CLIO bulk formulation
             ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
-            ! (zalb_ice) is computed within the bulk routine
-            CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, zalb_ice )
-            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            ! (alb_ice) is computed within the bulk routine
+                                 CALL blk_ice_clio_flx( t_su, zalb_cs, zalb_os, alb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE( jp_core )                                       ! CORE bulk formulation
             ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-            CALL blk_ice_core_flx( t_su, zalb_ice )
-            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL blk_ice_core_flx( t_su, alb_ice )
+            IF( ln_mixcpl      ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          CASE ( jp_purecpl )
             ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su )
-            ! clem: evap_ice is forced to 0 in coupled mode for now 
-            !       but it needs to be changed (along with modif in limthd_dh) once heat flux from evap will be avail. from atm. models
-            evap_ice  (:,:,:) = 0._wp   ;   devap_ice (:,:,:) = 0._wp
-            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+                                 CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+            IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          END SELECT
-         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs )
 
          !----------------------------!
@@ -264 +261 @@
       !!----------------------------------------------------------------------
       INTEGER :: ierr
+      INTEGER :: ji, jj
       !!----------------------------------------------------------------------
       IF(lwp) WRITE(numout,*)
@@ -320 +318 @@
       tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
       !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( gphit(ji,jj) > 0._wp ) THEN  ;  rn_amax_2d(ji,jj) = rn_amax_n  ! NH
+            ELSE                             ;  rn_amax_2d(ji,jj) = rn_amax_s  ! SH
+            ENDIF
+        ENDDO
+      ENDDO 
+      !
       nstart = numit  + nn_fsbc      
       nitrun = nitend - nit000 + 1 
@@ -342 +348 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir,  &
-         &                ln_limdyn, rn_amax, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
+         &                ln_limdyn, rn_amax_n, rn_amax_s, ln_limdiahsb, ln_limdiaout, ln_icectl, iiceprt, jiceprt  
       !!-------------------------------------------------------------------
       !                    
@@ -363 +369 @@
          WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
-         WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax 
+         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n 
+         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s
          WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
          WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
@@ -578 +585 @@
       sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-      sfx_res(:,:) = 0._wp
+      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
       
       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
@@ -594 +601 @@
       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp 
       hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
-      hfx_err_dif(:,:) = 0._wp   ;
-
+      hfx_err_dif(:,:) = 0._wp
+      wfx_err_sub(:,:) = 0._wp
+      
       afx_tot(:,:) = 0._wp   ;
       afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

@@ -150 +150 @@
 
          ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
-         tfu(:,:) = eos_fzp( sss_m ) +  rt0 
+         CALL eos_fzp( sss_m(:,:), tfu(:,:) )
+         tfu(:,:) = tfu(:,:) + rt0
 
          zsist (:,:,1) = sist (:,:) + rt0 * ( 1. - tmask(:,:,1) )

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

@@ -26 +26 @@
    USE zdfbfr
    USE fldread         ! read input field at current time step
-
-
+   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -53 +52 @@
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
-#if defined key_agrif
-   ! AGRIF can not handle these arrays as integers. The reason is a mystery but problems avoided by declaring them as reals
-   REAL(wp),    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-                                                                                          !: (first wet level and last level include in the tbl)
-#else
    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
-#endif
 
 
@@ -90 +83 @@
     INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
     INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
     REAL(wp)                     ::   rmin
     REAL(wp)                     ::   zhk
-    CHARACTER(len=256)           ::   cfisf, cvarzisf, cvarhisf   ! name for isf file
+    REAL(wp)                     ::   zt_frz, zpress
+    CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
     CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
     CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
     INTEGER           ::   ios           ! Local integer output status for namelist read
+
+    REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
+    REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
       !
       !!---------------------------------------------------------------------
@@ -176 +174 @@
               DO jj = 1, jpj
                   jk = 2
-                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
+                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. gdepw_0(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
                   misfkt(ji,jj) = jk-1
                END DO
@@ -194 +192 @@
          END IF
          
+         ! save initial top boundary layer thickness         
          rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
+
+      END IF
+
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
+         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
+         !
+      ENDIF
+
+      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
 
          ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
@@ -205 +216 @@
 
                ! determine the deepest level influenced by the boundary layer
-               ! test on tmask useless ?????
                DO jk = ikt, mbkt(ji,jj)
                   IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
@@ -217 +227 @@
             END DO
          END DO
-         
-      END IF
-
-      !                                            ! ---------------------------------------- !
-      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
-         !                                         ! ---------------------------------------- !
-         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
-         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
-         !
-      ENDIF
-
-      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
-
 
          ! compute salf and heat flux
@@ -256 +253 @@
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+
+            IF( lk_oasis) THEN
+            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
+
+              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
+              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
+              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
+ 
+              ! All related global sums must be done bit reproducibly
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+               ! check
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
+
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+      
+               ! check
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+            ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             stbl(:,:)   = soce
@@ -264 +302 @@
             !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
             fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+
+            IF( lk_oasis) THEN
+            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
+
+              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
+              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
+              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
+
+              ! All related global sums must be done bit reproducibly
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+               ! check
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
+
+               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+      
+               ! check
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
+
+               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+            ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
@@ -270 +349 @@
          END IF
          ! compute tsc due to isf
-         ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable).
-         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp !
+         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
+!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
+         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
          
          ! salt effect already take into account in vertical advection
          risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
-          
+
+         ! output
+         IF( iom_use('qlatisf' ) )   CALL iom_put('qlatisf', qisf)
+         IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce )
+
+         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
+         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
+ 
          ! lbclnk
          CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
@@ -281 +369 @@
          CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
          CALL lbc_lnk(qisf(:,:)     ,'T',1.)
+
+!=============================================================================================================================================
+         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
+            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
+            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
+
+            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
+            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
+            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
+            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
+
+            DO jj = 1,jpj
+               DO ji = 1,jpi
+                  ikt = misfkt(ji,jj)
+                  ikb = misfkb(ji,jj)
+                  DO jk = ikt, ikb - 1
+                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                  END DO
+                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
+                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
+                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
+               END DO
+            END DO
+
+            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
+            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
+            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
+            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
+
+            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
+            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
+         END IF
+!=============================================================================================================================================
 
          IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
@@ -295 +418 @@
          ENDIF
          ! 
-         ! output
-         CALL iom_put('qisf'  , qisf)
-         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
       END IF
   
@@ -370 +490 @@
              ! Calculate freezing temperature
                 zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
-                zt_frz = eos_fzp(tsb(ji,jj,ik,jp_sal), zpress) 
+                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
                 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
              ENDDO
@@ -452 +572 @@
       zti(:,:)=tinsitu( ttbl, stbl, zpress )
 ! Calculate freezing temperature
-      zfrz(:,:)=eos_fzp( sss_m(:,:), zpress )
+      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
 
       
@@ -472 +592 @@
 
                      nit = nit + 1
-                     IF (nit .GE. 100) THEN
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj
-                        !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx
-                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
-                     END IF
+                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+
 ! save gammat and compute zhtflx_b
                      zgammat2d(ji,jj)=zgammat
@@ -794 +911 @@
                ! test on tmask useless ?????
                DO jk = ikt, mbkt(ji,jj)
-!                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
                END DO
                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

@@ -179 +179 @@
 
       !                          ! Checks:
-      IF( nn_isf .EQ. 0 ) THEN                      ! no specific treatment in vicinity of ice shelf 
+      IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf 
          IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
-         fwfisf  (:,:) = 0.0_wp
-         fwfisf_b(:,:) = 0.0_wp
+         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
+         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
+         rdivisf       = 0.0_wp
       END IF
       IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero
@@ -265 +266 @@
       ENDIF
       !
-      IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
-      !                                                     !                                            (2) the use of nn_fsbc
+      IF( lk_oasis ) THEN
+         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )         
+         CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
+                                      !                                            (2) the use of nn_fsbc
+      ENDIF
 
 !     nn_fsbc initialization if OPA-SAS coupling via OASIS
@@ -339 +343 @@
          emp_b(:,:) = emp(:,:)
          sfx_b(:,:) = sfx(:,:)
+         IF ( ln_rnf ) THEN
+            rnf_b    (:,:  ) = rnf    (:,:  )
+            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+         ENDIF
       ENDIF
       !                                            ! ---------------------------------------- !
@@ -455 +463 @@
       !                                                ! ---------------------------------------- !
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         CALL iom_put( "empmr" , emp  - rnf )                   ! upward water flux
+         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux
+         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
          CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
                                                                 ! (includes virtual salt flux beneath ice 

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

@@ -52 +52 @@
    REAL(wp)                   ::   rn_hrnf         !: runoffs, depth over which enhanced vertical mixing is used
    REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s]
-   REAL(wp)                   ::   rn_rfact        !: multiplicative factor for runoff
+   REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff
 
    LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis
@@ -109 +109 @@
       !
       CALL wrk_alloc( jpi,jpj, ztfrz)
-
-      !                                            ! ---------------------------------------- !
-      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
-         !                                         ! ---------------------------------------- !
-         rnf_b    (:,:  ) = rnf    (:,:  )               ! Swap the ocean forcing fields except at nit000
-         rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)               ! where before fields are set at the end of the routine
-         !
-      ENDIF
-
+      !
       !                                            !-------------------!
       !                                            !   Update runoff   !
@@ -126 +118 @@
       IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
       !
-      ! Runoff reduction only associated to the ORCA2_LIM configuration
-      ! when reading the NetCDF file runoff_1m_nomask.nc
-      IF( cp_cfg == 'orca' .AND. jp_cfg == 2 .AND. .NOT. l_rnfcpl )   THEN
-         WHERE( 40._wp < gphit(:,:) .AND. gphit(:,:) < 65._wp )
-            sf_rnf(1)%fnow(:,:,1) = 0.85 * sf_rnf(1)%fnow(:,:,1)
-         END WHERE
-      ENDIF
-      !
       IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
          !
          IF( .NOT. l_rnfcpl )   rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt
+         CALL lbc_lnk(rnf(:,:), 'T', 1._wp)
          !
          !                                                     ! set temperature & salinity content of runoffs

branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90

@@ -31 +31 @@
 CONTAINS
 
-   SUBROUTINE upd_tide( kt, kit, kbaro, koffset )
+   SUBROUTINE upd_tide( kt, kit, time_offset )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE upd_tide  ***
@@ -42 +42 @@
       !!----------------------------------------------------------------------      
       INTEGER, INTENT(in)           ::   kt      ! ocean time-step index
-      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T only)
-      INTEGER, INTENT(in), OPTIONAL ::   kbaro   ! number of sub-time-step           (lk_dynspg_ts=T only)
-      INTEGER, INTENT(in), OPTIONAL ::   koffset ! time offset in number 
-                                                 ! of sub-time-steps                 (lk_dynspg_ts=T only)
+      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T)
+      INTEGER, INTENT(in), OPTIONAL ::   time_offset ! time offset in number 
+                                                     ! of internal steps             (lk_dynspg_ts=F)
+                                                     ! of external steps             (lk_dynspg_ts=T)
       !
       INTEGER  ::   joffset      ! local integer
@@ -57 +57 @@
       !
       joffset = 0
-      IF( PRESENT( koffset ) )   joffset = koffset
+      IF( PRESENT( time_offset ) )   joffset = time_offset
       !
-      IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   THEN
-         zt = zt + ( kit + 0.5_wp * ( joffset - 1 ) ) * rdt / REAL( kbaro, wp )
+      IF( PRESENT( kit ) )   THEN
+         zt = zt + ( kit +  joffset - 1 ) * rdt / REAL( nn_baro, wp )
       ELSE
          zt = zt + joffset * rdt
@@ -74 +74 @@
       IF( ln_tide_ramp ) THEN         ! linear increase if asked
          zt = ( kt - nit000 ) * rdt
-         IF( PRESENT( kit ) .AND. PRESENT( kbaro ) )   zt = zt + kit * rdt / REAL( kbaro, wp )
+         IF( PRESENT( kit ) )   zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp )
          zramp = MIN(  MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp  )
          pot_astro(:,:) = zramp * pot_astro(:,:)
@@ -86 +86 @@
   !!----------------------------------------------------------------------
 CONTAINS
-  SUBROUTINE upd_tide( kt, kit, kbaro, koffset )          ! Empty routine
+  SUBROUTINE upd_tide( kt, kit, time_offset )  ! Empty routine
     INTEGER, INTENT(in)           ::   kt      !  integer  arg, dummy routine
     INTEGER, INTENT(in), OPTIONAL ::   kit     !  optional arg, dummy routine
-    INTEGER, INTENT(in), OPTIONAL ::   kbaro   !  optional arg, dummy routine
-    INTEGER, INTENT(in), OPTIONAL ::   koffset !  optional arg, dummy routine
+    INTEGER, INTENT(in), OPTIONAL ::   time_offset !  optional arg, dummy routine
     WRITE(*,*) 'upd_tide: You should not have seen this print! error?', kt
   END SUBROUTINE upd_tide