Changeset 5407 for trunk/NEMOGCM/NEMO/OPA_SRC

Timestamp:

2015-06-11T21:13:22+02:00 (9 years ago)

Author:

smasson

Message:

merge dev_r5218_CNRS17_coupling into the trunk

Location:

trunk/NEMOGCM/NEMO/OPA_SRC

Files:

• DIA/diawri.F90 (modified) (6 diffs)
• DYN/divcur.F90 (modified) (1 diff)
• IOM/in_out_manager.F90 (modified) (1 diff)
• IOM/iom.F90 (modified) (7 diffs)
• IOM/restart.F90 (modified) (5 diffs)
• LBC/lib_mpp.F90 (modified) (4 diffs)
• SBC/cpl_oasis3.F90 (modified) (8 diffs)
• SBC/fldread.F90 (modified) (2 diffs)
• SBC/sbc_ice.F90 (modified) (5 diffs)
• SBC/sbc_oce.F90 (modified) (6 diffs)
• SBC/sbcblk_clio.F90 (modified) (16 diffs)
• SBC/sbcblk_core.F90 (modified) (11 diffs)
• SBC/sbccpl.F90 (modified) (47 diffs)
• SBC/sbcice_cice.F90 (modified) (8 diffs)
• SBC/sbcice_if.F90 (modified) (1 diff)
• SBC/sbcice_lim.F90 (modified) (15 diffs)
• SBC/sbcice_lim_2.F90 (modified) (8 diffs)
• SBC/sbcmod.F90 (modified) (17 diffs)
• SBC/sbcrnf.F90 (modified) (12 diffs)
• SBC/sbcssm.F90 (modified) (12 diffs)
• TRA/eosbn2.F90 (modified) (2 diffs)
• TRA/traqsr.F90 (modified) (12 diffs)
• ZDF/zdftke.F90 (modified) (1 diff)
• nemogcm.F90 (modified) (4 diffs)
• step.F90 (modified) (4 diffs)

trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

@@ -593 +593 @@
          ENDIF
 
-         IF( .NOT. lk_cpl ) THEN
+         IF( .NOT. ln_cpl ) THEN
             CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -602 +602 @@
          ENDIF
 
-         IF( lk_cpl .AND. nn_ice <= 1 ) THEN
+         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
             CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -625 +625 @@
 #endif
 
-         IF( lk_cpl .AND. nn_ice == 2 ) THEN
+         IF( ln_cpl .AND. nn_ice == 2 ) THEN
             CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -780 +780 @@
       ENDIF
 
-      IF( .NOT. lk_cpl ) THEN
+      IF( .NOT. ln_cpl ) THEN
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
@@ -786 +786 @@
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
-      IF( lk_cpl .AND. nn_ice <= 1 ) THEN
+      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
@@ -802 +802 @@
 #endif
 
-      IF( lk_cpl .AND. nn_ice == 2 ) THEN
+      IF( ln_cpl .AND. nn_ice == 2 ) THEN
          CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
          CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo

trunk/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

@@ -98 +98 @@
       !
       CALL wrk_alloc( jpi  , jpj+2, zwu               )
-      CALL wrk_alloc( jpi+4, jpj  , zwv, kjstart = -1 )
+      CALL wrk_alloc( jpi+4, jpj  , zwv, kistart = -1 )
       !
       IF( kt == nit000 ) THEN

trunk/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

@@ -149 +149 @@
    LOGICAL       ::   lwp      = .FALSE.    !: boolean : true on the 1st processor only .OR. ln_ctl
    LOGICAL       ::   lsp_area = .TRUE.     !: to make a control print over a specific area
+   CHARACTER(lc) ::   cxios_context         !: context name used in xios
 
    !!----------------------------------------------------------------------

trunk/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

@@ -121 +121 @@
       CALL set_scalar
 
-      IF( TRIM(cdname) == "nemo" ) THEN  
+      IF( TRIM(cdname) == TRIM(cxios_context) ) THEN  
          CALL set_grid( "T", glamt, gphit ) 
          CALL set_grid( "U", glamu, gphiu )
@@ -129 +129 @@
       ENDIF
 
-      IF( TRIM(cdname) == "nemo_crs" ) THEN  
+      IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN  
          CALL dom_grid_crs   ! Save the parent grid information  & Switch to coarse grid domain
          !
@@ -1212 +1212 @@
       CALL iom_swap( cdname )   ! swap to cdname context
       CALL xios_update_calendar(kt)
-      IF( cdname /= "nemo" ) CALL iom_swap( "nemo" )   ! return back to nemo context
+      IF( cdname /= TRIM(cxios_context) ) CALL iom_swap( TRIM(cxios_context) )   ! return back to nemo context
       !
    END SUBROUTINE iom_setkt
@@ -1222 +1222 @@
          CALL iom_swap( cdname )   ! swap to cdname context
          CALL xios_context_finalize() ! finalize the context
-         IF( cdname /= "nemo" ) CALL iom_swap( "nemo" )   ! return back to nemo context
+         IF( cdname /= TRIM(cxios_context) ) CALL iom_swap( TRIM(cxios_context) )   ! return back to nemo context
       ENDIF
       !
@@ -1291 +1291 @@
          CASE ('T', 'W')
             icnr = -1 ; jcnr = -1
-            IF( TRIM(cdname) == "nemo_crs" ) THEN
+            IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN
                z_cnr(:,:,1) = gphif_crs ; z_cnr(:,:,2) = glamf_crs
                z_pnt(:,:,1) = gphit_crs ; z_pnt(:,:,2) = glamt_crs
@@ -1300 +1300 @@
          CASE ('U')
             icnr =  0 ; jcnr = -1
-            IF( TRIM(cdname) == "nemo_crs" ) THEN
+            IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN
                z_cnr(:,:,1) = gphiv_crs ; z_cnr(:,:,2) = glamv_crs
                z_pnt(:,:,1) = gphiu_crs ; z_pnt(:,:,2) = glamu_crs
@@ -1309 +1309 @@
          CASE ('V')
             icnr = -1 ; jcnr =  0
-            IF( TRIM(cdname) == "nemo_crs" ) THEN
+            IF( TRIM(cdname) == TRIM(cxios_context)//"_crs" ) THEN
                z_cnr(:,:,1) = gphiu_crs ; z_cnr(:,:,2) = glamu_crs
                z_pnt(:,:,1) = gphiv_crs ; z_pnt(:,:,2) = glamv_crs

trunk/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

@@ -24 +24 @@
    USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
    USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE sbc_ice, ONLY : lk_lim3
 
    IMPLICIT NONE
@@ -135 +134 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      )
                      !
-      IF( lk_lim3 )  CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
-                     !
                      CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields
                      CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn        )
@@ -148 +145 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       )
 #endif
-                  IF( lk_lim3 ) THEN
-                     CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev'  , fraqsr_1lev     ) !clem modif
-                  ENDIF
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
@@ -236 +230 @@
          CALL iom_get( numror, jpdom_autoglo, 'hdivb'  , hdivb   )
          CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb    )
-         IF( lk_lim3 )   CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
       ELSE
          neuler = 0
@@ -279 +272 @@
          ENDIF
 
-         IF( lk_lim3 .AND. .NOT. lk_vvl ) THEN
-            DO jk = 1, jpk
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
-            END DO
-         ENDIF
-
-      ENDIF
-      !
-      IF( lk_lim3 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev )
       ENDIF
       !

trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

@@ -164 +164 @@
 
 
-   FUNCTION mynode( ldtxt, kumnam_ref , kumnam_cfg , kumond , kstop, localComm )
+   FUNCTION mynode( ldtxt, ldname, kumnam_ref , kumnam_cfg , kumond , kstop, localComm )
       !!----------------------------------------------------------------------
       !!                  ***  routine mynode  ***
@@ -171 +171 @@
       !!----------------------------------------------------------------------
       CHARACTER(len=*),DIMENSION(:), INTENT(  out) ::   ldtxt
+      CHARACTER(len=*)             , INTENT(in   ) ::   ldname
       INTEGER                      , INTENT(in   ) ::   kumnam_ref     ! logical unit for reference namelist
       INTEGER                      , INTENT(in   ) ::   kumnam_cfg     ! logical unit for configuration namelist
@@ -297 +298 @@
 
       IF( mynode == 0 ) THEN
-        CALL ctl_opn( kumond, 'output.namelist.dyn', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. , 1 )
-        WRITE(kumond, nammpp)      
+         CALL ctl_opn( kumond, TRIM(ldname), 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. , 1 )
+         WRITE(kumond, nammpp)      
       ENDIF
       !
@@ -3192 +3193 @@
    END FUNCTION lib_mpp_alloc
 
-   FUNCTION mynode( ldtxt, kumnam_ref, knumnam_cfg,  kumond , kstop, localComm ) RESULT (function_value)
+   FUNCTION mynode( ldtxt, ldname, kumnam_ref, knumnam_cfg,  kumond , kstop, localComm ) RESULT (function_value)
       INTEGER, OPTIONAL            , INTENT(in   ) ::   localComm
       CHARACTER(len=*),DIMENSION(:) ::   ldtxt
+      CHARACTER(len=*) ::   ldname
       INTEGER ::   kumnam_ref, knumnam_cfg , kumond , kstop
       IF( PRESENT( localComm ) .OR. .NOT.PRESENT( localComm ) )   function_value = 0
       IF( .FALSE. )   ldtxt(:) = 'never done'
-      CALL ctl_opn( kumond, 'output.namelist.dyn', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. , 1 )
+      CALL ctl_opn( kumond, TRIM(ldname), 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. , 1 )
    END FUNCTION mynode
 

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90

@@ -15 +15 @@
    !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
-   !!   'key_oasis3'                    coupled Ocean/Atmosphere via OASIS3
+   !!   'key_oasis3'                    coupled Ocean/Atmosphere via OASIS3-MCT
+   !!   'key_oa3mct_v3'                 to be added for OASIS3-MCT version 3
    !!----------------------------------------------------------------------
    !!   cpl_init     : initialization of coupled mode communication
@@ -61 +62 @@
 #endif
 
-   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=40        ! Maximum number of coupling fields
+   INTEGER                    ::   nrcv         ! total number of fields received 
+   INTEGER                    ::   nsnd         ! total number of fields sent 
+   INTEGER                    ::   ncplmodel    ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
+   INTEGER, PUBLIC, PARAMETER ::   nmaxfld=50   ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields
    INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields
@@ -86 +90 @@
 CONTAINS
 
-   SUBROUTINE cpl_init( kl_comm )
+   SUBROUTINE cpl_init( cd_modname, kl_comm )
       !!-------------------------------------------------------------------
       !!             ***  ROUTINE cpl_init  ***
@@ -95 +99 @@
       !! ** Method  :   OASIS3 MPI communication 
       !!--------------------------------------------------------------------
-      INTEGER, INTENT(out) ::   kl_comm   ! local communicator of the model
+      CHARACTER(len = *), INTENT(in) ::   cd_modname   ! model name as set in namcouple file
+      INTEGER          , INTENT(out) ::   kl_comm      ! local communicator of the model
       !!--------------------------------------------------------------------
 
@@ -104 +109 @@
       ! 1st Initialize the OASIS system for the application
       !------------------------------------------------------------------
-      CALL oasis_init_comp ( ncomp_id, 'oceanx', nerror )
+      CALL oasis_init_comp ( ncomp_id, TRIM(cd_modname), nerror )
       IF ( nerror /= OASIS_Ok ) &
          CALL oasis_abort (ncomp_id, 'cpl_init', 'Failure in oasis_init_comp')
@@ -144 +149 @@
       IF(lwp) WRITE(numout,*)
 
+      ncplmodel = kcplmodel
       IF( kcplmodel > nmaxcpl ) THEN
-         CALL oasis_abort ( ncomp_id, 'cpl_define', 'kcplmodel is larger than nmaxcpl, increase nmaxcpl')   ;   RETURN
+         CALL oasis_abort ( ncomp_id, 'cpl_define', 'ncplmodel is larger than nmaxcpl, increase nmaxcpl')   ;   RETURN
       ENDIF
+
+      nrcv = krcv
+      IF( nrcv > nmaxfld ) THEN
+         CALL oasis_abort ( ncomp_id, 'cpl_define', 'nrcv is larger than nmaxfld, increase nmaxfld')   ;   RETURN
+      ENDIF
+
+      nsnd = ksnd
+      IF( nsnd > nmaxfld ) THEN
+         CALL oasis_abort ( ncomp_id, 'cpl_define', 'nsnd is larger than nmaxfld, increase nmaxfld')   ;   RETURN
+      ENDIF
+
       !
       ! ... Define the shape for the area that excludes the halo
@@ -400 +417 @@
 
 
-   INTEGER FUNCTION cpl_freq( kid )  
+   INTEGER FUNCTION cpl_freq( cdfieldname )  
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_freq  ***
@@ -406 +423 @@
       !! ** Purpose : - send back the coupling frequency for a particular field
       !!----------------------------------------------------------------------
-      INTEGER,INTENT(in) ::   kid   ! variable index
-      !!
+      CHARACTER(len = *), INTENT(in) ::   cdfieldname    ! field name as set in namcouple file
+      !!
+      INTEGER               :: id
       INTEGER               :: info
       INTEGER, DIMENSION(1) :: itmp
+      INTEGER               :: ji,jm     ! local loop index
+      INTEGER               :: mop
       !!----------------------------------------------------------------------
-      CALL oasis_get_freqs(kid, 1, itmp, info)
-      cpl_freq = itmp(1)
+      cpl_freq = 0   ! defaut definition
+      id = -1        ! defaut definition
+      !
+      DO ji = 1, nsnd
+         IF (ssnd(ji)%laction ) THEN
+            DO jm = 1, ncplmodel
+               IF( ssnd(ji)%nid(1,jm) /= -1 ) THEN
+                  IF( TRIM(cdfieldname) == TRIM(ssnd(ji)%clname) ) THEN
+                     id = ssnd(ji)%nid(1,jm)
+                     mop = OASIS_Out
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDIF
+      ENDDO
+      DO ji = 1, nrcv
+         IF (srcv(ji)%laction ) THEN
+            DO jm = 1, ncplmodel
+               IF( srcv(ji)%nid(1,jm) /= -1 ) THEN
+                  IF( TRIM(cdfieldname) == TRIM(srcv(ji)%clname) ) THEN
+                     id = srcv(ji)%nid(1,jm)
+                     mop = OASIS_In
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDIF
+      ENDDO
+      !
+      IF( id /= -1 ) THEN
+#if defined key_oa3mct_v3
+         CALL oasis_get_freqs(id, mop, 1, itmp, info)
+#else
+         CALL oasis_get_freqs(id,      1, itmp, info)
+#endif
+         cpl_freq = itmp(1)
+      ENDIF
       !
    END FUNCTION cpl_freq

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

@@ -154 +154 @@
       IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1
 
-      it_offset = 0
+      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
+      ELSE                                      ;   it_offset = 0
+      ENDIF
       IF( PRESENT(kt_offset) )   it_offset = kt_offset
 
@@ -452 +454 @@
       ENDIF
       !
-      it_offset = 0
+      IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc
+      ELSE                                      ;   it_offset = 0
+      ENDIF
       IF( PRESENT(kt_offset) )   it_offset = kt_offset
       IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) )

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

@@ -68 +68 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice            [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt            !: category topmelt
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   botmelt            !: category botmelt
+
+#if defined  key_lim3
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   evap_ice       !: sublimation                              [kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   devap_ice      !: sublimation sensitivity                [kg/m2/s/K]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qns_oce        !: non solar heat flux over ocean              [W/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qsr_oce        !: non solar heat flux over ocean              [W/m2]
+   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(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s]
+#endif
+#if defined key_lim3 || defined key_lim2
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s]
+#endif
 
 #if defined key_cice
@@ -99 +113 @@
 #endif
 
-#if defined key_lim3 || defined key_cice
-   ! not used with LIM2
+#if defined key_cice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tatm_ice       !: air temperature [K]
 #endif
@@ -124 +137 @@
       ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
          &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
-         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) ,     &
-         &      alb_ice (jpi,jpj,jpl) ,                             &
-         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     ,     &
+         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   &
+         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   &
          &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     &
-#if defined key_lim3
-         &      tatm_ice(jpi,jpj)     ,                             &
-#endif
 #if defined key_lim2
          &      a_i(jpi,jpj,jpl)      ,                             &
+#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)  ,  &
 #endif
          &      emp_ice(jpi,jpj)      ,  STAT= ierr(1) )
@@ -144 +158 @@
                 a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , &
                 STAT= ierr(1) )
-      IF( lk_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , &
+      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , &
          &                     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)   , &
@@ -152 +166 @@
          !
 #if defined key_cice || defined key_lim2
-      IF( lk_cpl )   ALLOCATE( ht_i(jpi,jpj,jpl) , ht_s(jpi,jpj,jpl) , STAT=ierr(5) )
+      IF( ln_cpl )   ALLOCATE( ht_i(jpi,jpj,jpl) , ht_s(jpi,jpj,jpl) , STAT=ierr(5) )
 #endif
 

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

@@ -36 +36 @@
    LOGICAL , PUBLIC ::   ln_blk_mfs     !: MFS  bulk formulation
 #if defined key_oasis3
-   LOGICAL , PUBLIC ::   lk_cpl = .TRUE.  !: coupled formulation
+   LOGICAL , PUBLIC ::   lk_oasis = .TRUE.  !: OASIS used
 #else
-   LOGICAL , PUBLIC ::   lk_cpl = .FALSE. !: coupled formulation
-#endif
+   LOGICAL , PUBLIC ::   lk_oasis = .FALSE. !: OASIS unused
+#endif
+   LOGICAL , PUBLIC ::   ln_cpl         !: ocean-atmosphere coupled formulation
+   LOGICAL , PUBLIC ::   ln_mixcpl      !: ocean-atmosphere forced-coupled mixed formulation
    LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr)
    LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths
@@ -50 +52 @@
    !                                             !: =1 levitating ice with mass and salt exchange but no presure effect
    !                                             !: =2 embedded sea-ice (full salt and mass exchanges and pressure)
-   INTEGER , PUBLIC :: nn_limflx        !: LIM3 Multi-category heat flux formulation
+   INTEGER , PUBLIC ::   nn_components  !: flag for sbc module (including sea-ice) coupling mode (see component definition below) 
+   INTEGER , PUBLIC ::   nn_limflx      !: LIM3 Multi-category heat flux formulation
    !                                             !: =-1  Use of per-category fluxes
    !                                             !: = 0  Average per-category fluxes
@@ -69 +72 @@
    !!           switch definition (improve readability)
    !!----------------------------------------------------------------------
-   INTEGER , PUBLIC, PARAMETER ::   jp_gyre    = 0        !: GYRE analytical formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_ana     = 1        !: analytical      formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_flx     = 2        !: flux            formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_clio    = 3        !: CLIO bulk       formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_core    = 4        !: CORE bulk       formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_cpl     = 5        !: Coupled         formulation
-   INTEGER , PUBLIC, PARAMETER ::   jp_mfs     = 6        !: MFS  bulk       formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_gyre    = 0        !: GYRE analytical               formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_ana     = 1        !: analytical                    formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_flx     = 2        !: flux                          formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_clio    = 3        !: CLIO bulk                     formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_core    = 4        !: CORE bulk                     formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_purecpl = 5        !: Pure ocean-atmosphere Coupled formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_mfs     = 6        !: MFS  bulk                     formulation
+   INTEGER , PUBLIC, PARAMETER ::   jp_none    = 7        !: for OPA when doing coupling via SAS module
    INTEGER , PUBLIC, PARAMETER ::   jp_esopa   = -1       !: esopa test, ALL formulations
    
+   !!----------------------------------------------------------------------
+   !!           component definition
+   !!----------------------------------------------------------------------
+   INTEGER , PUBLIC, PARAMETER ::   jp_iam_nemo = 0      !: Initial single executable configuration 
+                                                         !  (no internal OASIS coupling)
+   INTEGER , PUBLIC, PARAMETER ::   jp_iam_opa  = 1      !: Multi executable configuration - OPA component
+                                                         !  (internal OASIS coupling)
+   INTEGER , PUBLIC, PARAMETER ::   jp_iam_sas  = 2      !: Multi executable configuration - SAS component
+                                                         !  (internal OASIS coupling)
    !!----------------------------------------------------------------------
    !!              Ocean Surface Boundary Condition fields
@@ -111 +124 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm]
 #endif
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
 
    !!----------------------------------------------------------------------
@@ -122 +136 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m]
    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 [-]
 
    !! * Substitutions
@@ -155 +170 @@
          &      atm_co2(jpi,jpj) ,                                        &
 #endif
-         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       &
-         &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
+         &      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) )
          !
 #if defined key_vvl

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90

@@ -34 +34 @@
    USE albedo
    USE prtctl          ! Print control
-#if defined key_lim3
+#if defined key_lim3 
    USE ice
    USE sbc_ice         ! Surface boundary condition: ice fields
+   USE limthd_dh       ! for CALL lim_thd_snwblow
 #elif defined key_lim2
    USE ice_2
+   USE sbc_ice         ! Surface boundary condition: ice fields
+   USE par_ice_2       ! Surface boundary condition: ice fields
 #endif
 
@@ -45 +48 @@
 
    PUBLIC sbc_blk_clio        ! routine called by sbcmod.F90 
-   PUBLIC blk_ice_clio        ! routine called by sbcice_lim.F90 
+#if defined key_lim2 || defined key_lim3
+   PUBLIC blk_ice_clio_tau    ! routine called by sbcice_lim.F90 
+   PUBLIC blk_ice_clio_flx    ! routine called by sbcice_lim.F90 
+#endif
 
    INTEGER , PARAMETER ::   jpfld   = 7           ! maximum number of files to read 
@@ -378 +384 @@
          &     + sf(jp_prec)%fnow(:,:,1) * sf(jp_tair)%fnow(:,:,1) * zcprec   ! add    precip. heat content at Tair in Celcius
       qns(:,:) = qns(:,:) * tmask(:,:,1)
+#if defined key_lim3
+      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)
+      qsr_oce(:,:) = qsr(:,:)
+#endif
       ! NB: if sea-ice model, the snow precip are computed and the associated heat is added to qns (see blk_ice_clio)
 
-      CALL iom_put( "qlw_oce",   zqlw )   ! output downward longwave  heat over the ocean
-      CALL iom_put( "qsb_oce", - zqsb )   ! output downward sensible  heat over the ocean
-      CALL iom_put( "qla_oce", - zqla )   ! output downward latent    heat over the ocean
-      CALL iom_put( "qns_oce",   qns  )   ! output downward non solar heat over the ocean
+      IF ( nn_ice == 0 ) THEN
+         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave  heat over the ocean
+         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible  heat over the ocean
+         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent    heat over the ocean
+         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
+         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
+         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
+      ENDIF
 
       IF(ln_ctl) THEN
@@ -399 +414 @@
    END SUBROUTINE blk_oce_clio
 
-
-   SUBROUTINE blk_ice_clio(  pst   , palb_cs, palb_os, palb,  &
-      &                      p_taui, p_tauj, p_qns , p_qsr,   &
-      &                      p_qla , p_dqns, p_dqla,          &
-      &                      p_tpr , p_spr ,                  &
-      &                      p_fr1 , p_fr2 , cd_grid, pdim  )
+# if defined key_lim2 || defined key_lim3
+   SUBROUTINE blk_ice_clio_tau
       !!---------------------------------------------------------------------------
-      !!                     ***  ROUTINE blk_ice_clio  ***
+      !!                     ***  ROUTINE blk_ice_clio_tau  ***
+      !!                 
+      !!  ** Purpose :   Computation momentum flux at the ice-atm interface  
+      !!         
+      !!  ** Method  :   Read utau from a forcing file. Rearrange if C-grid
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   zcoef
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio_tau')
+
+      SELECT CASE( cp_ice_msh )
+
+      CASE( 'C' )                          ! C-grid ice dynamics
+
+         zcoef  = cai / cao                         ! Change from air-sea stress to air-ice stress
+         utau_ice(:,:) = zcoef * utau(:,:)
+         vtau_ice(:,:) = zcoef * vtau(:,:)
+
+      CASE( 'I' )                          ! I-grid ice dynamics:  I-point (i.e. F-point lower-left corner)
+
+         zcoef  = 0.5_wp * cai / cao                ! Change from air-sea stress to air-ice stress
+         DO jj = 2, jpj         ! stress from ocean U- and V-points to ice U,V point
+            DO ji = 2, jpi   ! I-grid : no vector opt.
+               utau_ice(ji,jj) = zcoef * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) )
+               vtau_ice(ji,jj) = zcoef * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) )
+            END DO
+         END DO
+
+         CALL lbc_lnk( utau_ice(:,:), 'I', -1. )   ;   CALL lbc_lnk( vtau_ice(:,:), 'I', -1. )   ! I-point
+
+      END SELECT
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice_clio: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ')
+      ENDIF
+
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio_tau')
+
+   END SUBROUTINE blk_ice_clio_tau
+#endif
+
+# if defined key_lim2 || defined key_lim3
+   SUBROUTINE blk_ice_clio_flx(  ptsu , palb_cs, palb_os, palb )
+      !!---------------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_clio_flx ***
       !!                 
       !!  ** Purpose :   Computation of the heat fluxes at ocean and snow/ice
@@ -428 +486 @@
       !!                         to take into account solid precip latent heat flux
       !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature                   [Kelvin]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   ptsu      ! ice surface temperature                   [Kelvin]
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_cs  ! ice albedo (clear    sky) (alb_ice_cs)         [-]
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb_os  ! ice albedo (overcast sky) (alb_ice_os)         [-]
       REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   palb     ! ice albedo (actual value)                      [-]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_taui   ! surface ice stress at I-point (i-component) [N/m2]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tauj   ! surface ice stress at I-point (j-component) [N/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qns    ! non solar heat flux over ice (T-point)      [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qsr    !     solar heat flux over ice (T-point)      [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qla    ! latent    heat flux over ice (T-point)      [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqns   ! non solar heat sensistivity  (T-point)      [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqla   ! latent    heat sensistivity  (T-point)      [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tpr    ! total precipitation          (T-point)   [Kg/m2/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_spr    ! solid precipitation          (T-point)   [Kg/m2/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr1    ! 1sr fraction of qsr penetration in ice         [-]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr2    ! 2nd fraction of qsr penetration in ice         [-]
-      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! type of sea-ice grid ("C" or "B" grid)
-      INTEGER, INTENT(in   )                      ::   pdim     ! number of ice categories
       !!
       INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
-      !!
-      REAL(wp) ::   zcoef, zmt1, zmt2, zmt3, ztatm3     ! temporary scalars
+      !!
+      REAL(wp) ::   zmt1, zmt2, zmt3, ztatm3                    ! temporary scalars
       REAL(wp) ::   ztaevbk, zind1, zind2, zind3, ztamr         !    -         -
       REAL(wp) ::   zesi, zqsati, zdesidt                       !    -         -
@@ -455 +499 @@
       REAL(wp) ::   zcshi, zclei, zrhovaclei, zrhovacshi        !    -         -
       REAL(wp) ::   ztice3, zticemb, zticemb2, zdqlw, zdqsb     !    -         -
+      REAL(wp) ::   z1_lsub                                     !    -         -
       !!
       REAL(wp), DIMENSION(:,:)  , POINTER ::   ztatm   ! Tair in Kelvin
@@ -461 +506 @@
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa   ! air density
       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw, z_qsb
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw
       !!---------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio')
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_clio_flx')
       !
       CALL wrk_alloc( jpi,jpj, ztatm, zqatm, zevsqr, zrhoa )
-      CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb )
-
-      ijpl  = pdim                           ! number of ice categories
+      CALL wrk_alloc( jpi,jpj, jpl, z_qlw, z_qsb )
+
       zpatm = 101000.                        ! atmospheric pressure  (assumed constant  here)
-
-#if defined key_lim3      
-      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init
-#endif
-      !                                                        ! surface ocean fluxes computed with CLIO bulk formulea
-      !------------------------------------!
-      !   momentum fluxes  (utau, vtau )   !
-      !------------------------------------!
-
-      SELECT CASE( cd_grid )
-      CASE( 'C' )                          ! C-grid ice dynamics
-         zcoef  = cai / cao                         ! Change from air-sea stress to air-ice stress
-         p_taui(:,:) = zcoef * utau(:,:)
-         p_tauj(:,:) = zcoef * vtau(:,:)
-      CASE( 'I' )                          ! I-grid ice dynamics:  I-point (i.e. F-point lower-left corner)
-         zcoef  = 0.5_wp * cai / cao                ! Change from air-sea stress to air-ice stress
-         DO jj = 2, jpj         ! stress from ocean U- and V-points to ice U,V point
-            DO ji = 2, jpi   ! I-grid : no vector opt.
-               p_taui(ji,jj) = zcoef * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) )
-               p_tauj(ji,jj) = zcoef * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) )
-            END DO
-         END DO
-         CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ;   CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point
-      END SELECT
-
-
+      !--------------------------------------------------------------------------------
       !  Determine cloud optical depths as a function of latitude (Chou et al., 1981).
       !  and the correction factor for taking into account  the effect of clouds 
-      !------------------------------------------------------
+      !--------------------------------------------------------------------------------
+
 !CDIR NOVERRCHK
 !CDIR COLLAPSE
@@ -525 +546 @@
             zmt2  = ( 272.0 - ztatm(ji,jj) ) / 38.0   ;   zind2 = MAX( 0.e0, SIGN( 1.e0, zmt2 ) )
             zmt3  = ( 281.0 - ztatm(ji,jj) ) / 18.0   ;   zind3 = MAX( 0.e0, SIGN( 1.e0, zmt3 ) )
-            p_spr(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) / rday   &      ! rday = converte mm/day to kg/m2/s
+            sprecip(ji,jj) = sf(jp_prec)%fnow(ji,jj,1) / rday   &      ! rday = converte mm/day to kg/m2/s
                &         * (          zind1      &                   ! solid  (snow) precipitation [kg/m2/s]
                &            + ( 1.0 - zind1 ) * (          zind2   * ( 0.5 + zmt2 )   &
@@ -535 +556 @@
             ! fraction of qsr_ice which is NOT absorbed in the thin surface layer
             ! and thus which penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
-            p_fr1(ji,jj) = 0.18  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj,1) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj,1) 
-            p_fr2(ji,jj) = 0.82  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj,1) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj,1)
-         END DO
-      END DO
-      CALL iom_put( 'snowpre', p_spr )   ! Snow precipitation 
+            fr1_i0(ji,jj) = 0.18  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj,1) ) + 0.35 * sf(jp_ccov)%fnow(ji,jj,1) 
+            fr2_i0(ji,jj) = 0.82  * ( 1.e0 - sf(jp_ccov)%fnow(ji,jj,1) ) + 0.65 * sf(jp_ccov)%fnow(ji,jj,1)
+         END DO
+      END DO
+      CALL iom_put( 'snowpre', sprecip )   ! Snow precipitation 
       
       !-----------------------------------------------------------!
       !  snow/ice Shortwave radiation   (abedo already computed)  !
       !-----------------------------------------------------------!
-      CALL blk_clio_qsr_ice( palb_cs, palb_os, p_qsr )
-      
-      DO jl = 1, ijpl
+      CALL blk_clio_qsr_ice( palb_cs, palb_os, qsr_ice )
+      
+      DO jl = 1, jpl
          palb(:,:,jl) = ( palb_cs(:,:,jl) * ( 1.e0 - sf(jp_ccov)%fnow(:,:,1) )   &
             &         +   palb_os(:,:,jl) * sf(jp_ccov)%fnow(:,:,1) )
@@ -552 +573 @@
 
       !                                     ! ========================== !
-      DO jl = 1, ijpl                       !  Loop over ice categories  !
+      DO jl = 1, jpl                       !  Loop over ice categories  !
          !                                  ! ========================== !
 !CDIR NOVERRCHK
@@ -566 +587 @@
                ztaevbk = ztatm3 * ztatm(ji,jj) * zcldeff * ( 0.39 - 0.05 * zevsqr(ji,jj) ) 
                !
-               z_qlw(ji,jj,jl) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( pst(ji,jj,jl) - ztatm(ji,jj) ) ) 
+               z_qlw(ji,jj,jl) = - emic * stefan * ( ztaevbk + 4. * ztatm3 * ( ptsu(ji,jj,jl) - ztatm(ji,jj) ) ) 
 
                !----------------------------------------
@@ -573 +594 @@
 
                ! vapour pressure at saturation of ice (tmask to avoid overflow in the exponential)
-               zesi =  611.0 * EXP( 21.8745587 * tmask(ji,jj,1) * ( pst(ji,jj,jl) - rtt )/ ( pst(ji,jj,jl) - 7.66 ) )
+               zesi =  611.0 * EXP( 21.8745587 * tmask(ji,jj,1) * ( ptsu(ji,jj,jl) - rtt )/ ( ptsu(ji,jj,jl) - 7.66 ) )
                ! humidity close to the ice surface (at saturation)
                zqsati   = ( 0.622 * zesi ) / ( zpatm - 0.378 * zesi )
                
                !  computation of intermediate values
-               zticemb  = pst(ji,jj,jl) - 7.66
+               zticemb  = ptsu(ji,jj,jl) - 7.66
                zticemb2 = zticemb * zticemb  
-               ztice3   = pst(ji,jj,jl) * pst(ji,jj,jl) * pst(ji,jj,jl)
+               ztice3   = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
                zdesidt  = zesi * ( 9.5 * LOG( 10.0 ) * ( rtt - 7.66 )  / zticemb2 )
                
@@ -593 +614 @@
             
                !  sensible heat flux
-               z_qsb(ji,jj,jl) = zrhovacshi * ( pst(ji,jj,jl) - ztatm(ji,jj) )
+               z_qsb(ji,jj,jl) = zrhovacshi * ( ptsu(ji,jj,jl) - ztatm(ji,jj) )
             
                !  latent heat flux 
-               p_qla(ji,jj,jl) = MAX(  0.e0, zrhovaclei * ( zqsati - zqatm(ji,jj) )  )
+               qla_ice(ji,jj,jl) = MAX(  0.e0, zrhovaclei * ( zqsati - zqatm(ji,jj) )  )
               
                !  sensitivity of non solar fluxes (dQ/dT) (long-wave, sensible and latent fluxes)
@@ -603 +624 @@
                zdqla = zrhovaclei * ( zdesidt * ( zqsati * zqsati / ( zesi * zesi ) ) * ( zpatm / 0.622 ) )   
                !
-               p_dqla(ji,jj,jl) = zdqla                           ! latent flux sensitivity
-               p_dqns(ji,jj,jl) = -( zdqlw + zdqsb + zdqla )      !  total non solar sensitivity
+               dqla_ice(ji,jj,jl) = zdqla                           ! latent flux sensitivity
+               dqns_ice(ji,jj,jl) = -( zdqlw + zdqsb + zdqla )      !  total non solar sensitivity
             END DO
             !
@@ -616 +637 @@
       !
 !CDIR COLLAPSE
-      p_qns(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - p_qla (:,:,:)      ! Downward Non Solar flux
-!CDIR COLLAPSE
-      p_tpr(:,:)   = sf(jp_prec)%fnow(:,:,1) / rday                     ! total precipitation [kg/m2/s]
+      qns_ice(:,:,:) = z_qlw (:,:,:) - z_qsb (:,:,:) - qla_ice (:,:,:)      ! Downward Non Solar flux
+!CDIR COLLAPSE
+      tprecip(:,:)   = sf(jp_prec)%fnow(:,:,1) / rday                     ! total precipitation [kg/m2/s]
       !
       ! ----------------------------------------------------------------------------- !
@@ -625 +646 @@
 !CDIR COLLAPSE
       qns(:,:) = qns(:,:)                                                           &   ! update the non-solar heat flux with:
-         &     - p_spr(:,:) * lfus                                                  &   ! remove melting solid precip
-         &     + p_spr(:,:) * MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow - rt0 ) * cpic &   ! add solid P at least below melting
-         &     - p_spr(:,:) * sf(jp_tair)%fnow(:,:,1)                        * rcp      ! remove solid precip. at Tair
-      !
+         &     - sprecip(:,:) * lfus                                                  &   ! remove melting solid precip
+         &     + sprecip(:,:) * MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow - rt0 ) * cpic &   ! add solid P at least below melting
+         &     - sprecip(:,:) * sf(jp_tair)%fnow(:,:,1)                        * rcp      ! remove solid precip. at Tair
+
+#if defined key_lim3
+      ! ----------------------------------------------------------------------------- !
+      !    Distribute evapo, precip & associated heat over ice and ocean
+      ! ---------------=====--------------------------------------------------------- !
+      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
+
+      ! --- 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
+
+      ! --- evaporation minus precipitation --- !
+      CALL lim_thd_snwblow( pfrld, zsnw )          ! snow redistribution by wind
+      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * ( 1._wp - zsnw )
+      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
+      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
+
+      ! --- heat flux associated with emp --- !
+      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap
+         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip
+         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+
+      ! --- total solar and non solar fluxes --- !
+      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
+      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
+
+      ! --- 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 )
+
+      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
+#endif
+
 !!gm : not necessary as all input data are lbc_lnk...
-      CALL lbc_lnk( p_fr1  (:,:) , 'T', 1. )
-      CALL lbc_lnk( p_fr2  (:,:) , 'T', 1. )
-      DO jl = 1, ijpl
-         CALL lbc_lnk( p_qns (:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( p_dqns(:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( p_qla (:,:,jl) , 'T', 1. )
-         CALL lbc_lnk( p_dqla(:,:,jl) , 'T', 1. )
+      CALL lbc_lnk( fr1_i0  (:,:) , 'T', 1. )
+      CALL lbc_lnk( fr2_i0  (:,:) , 'T', 1. )
+      DO jl = 1, jpl
+         CALL lbc_lnk( qns_ice (:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( dqns_ice(:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( qla_ice (:,:,jl) , 'T', 1. )
+         CALL lbc_lnk( dqla_ice(:,:,jl) , 'T', 1. )
       END DO
 
 !!gm : mask is not required on forcing
-      DO jl = 1, ijpl
-         p_qns (:,:,jl) = p_qns (:,:,jl) * tmask(:,:,1)
-         p_qla (:,:,jl) = p_qla (:,:,jl) * tmask(:,:,1)
-         p_dqns(:,:,jl) = p_dqns(:,:,jl) * tmask(:,:,1)
-         p_dqla(:,:,jl) = p_dqla(:,:,jl) * tmask(:,:,1)
-      END DO
+      DO jl = 1, jpl
+         qns_ice (:,:,jl) = qns_ice (:,:,jl) * tmask(:,:,1)
+         qla_ice (:,:,jl) = qla_ice (:,:,jl) * tmask(:,:,1)
+         dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * tmask(:,:,1)
+         dqla_ice(:,:,jl) = dqla_ice(:,:,jl) * tmask(:,:,1)
+      END DO
+
+      CALL wrk_dealloc( jpi,jpj, ztatm, zqatm, zevsqr, zrhoa )
+      CALL wrk_dealloc( jpi,jpj, jpl  , z_qlw, z_qsb )
 
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab3d_1=z_qsb  , clinfo1=' blk_ice_clio: z_qsb  : ', tab3d_2=z_qlw  , clinfo2=' z_qlw  : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=p_qla  , clinfo1=' blk_ice_clio: z_qla  : ', tab3d_2=p_qsr  , clinfo2=' p_qsr  : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=p_dqns , clinfo1=' blk_ice_clio: p_dqns : ', tab3d_2=p_qns  , clinfo2=' p_qns  : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=p_dqla , clinfo1=' blk_ice_clio: p_dqla : ', tab3d_2=pst    , clinfo2=' pst    : ', kdim=ijpl)
-         CALL prt_ctl(tab2d_1=p_tpr  , clinfo1=' blk_ice_clio: p_tpr  : ', tab2d_2=p_spr  , clinfo2=' p_spr  : ')
-         CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_clio: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ')
+         CALL prt_ctl(tab3d_1=z_qsb  , clinfo1=' blk_ice_clio: z_qsb  : ', tab3d_2=z_qlw  , clinfo2=' z_qlw  : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=qla_ice  , clinfo1=' blk_ice_clio: z_qla  : ', tab3d_2=qsr_ice  , clinfo2=' qsr_ice  : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=dqns_ice , clinfo1=' blk_ice_clio: dqns_ice : ', tab3d_2=qns_ice  , clinfo2=' qns_ice  : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=dqla_ice , clinfo1=' blk_ice_clio: dqla_ice : ', tab3d_2=ptsu    , clinfo2=' ptsu    : ', kdim=jpl)
+         CALL prt_ctl(tab2d_1=tprecip  , clinfo1=' blk_ice_clio: tprecip  : ', tab2d_2=sprecip  , clinfo2=' sprecip  : ')
       ENDIF
 
-      CALL wrk_dealloc( jpi,jpj, ztatm, zqatm, zevsqr, zrhoa )
-      CALL wrk_dealloc( jpi,jpj,pdim, z_qlw, z_qsb )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio')
-      !
-   END SUBROUTINE blk_ice_clio
-
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_clio_flx')
+      !
+   END SUBROUTINE blk_ice_clio_flx
+
+#endif
 
    SUBROUTINE blk_clio_qsr_oce( pqsr_oce )

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -44 +44 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE lib_fortran     ! to use key_nosignedzero
+#if defined key_lim3
+   USE ice, ONLY       : u_ice, v_ice, jpl, pfrld, a_i_b
+   USE limthd_dh       ! for CALL lim_thd_snwblow
+#elif defined key_lim2
+   USE ice_2, ONLY     : u_ice, v_ice
+   USE par_ice_2
+#endif
 
    IMPLICIT NONE
@@ -49 +56 @@
 
    PUBLIC   sbc_blk_core         ! routine called in sbcmod module
-   PUBLIC   blk_ice_core         ! routine called in sbc_ice_lim module
+#if defined key_lim2 || defined key_lim3
+   PUBLIC   blk_ice_core_tau     ! routine called in sbc_ice_lim module
+   PUBLIC   blk_ice_core_flx     ! routine called in sbc_ice_lim module
+#endif
    PUBLIC   turb_core_2z         ! routine calles in sbcblk_mfs module
 
@@ -371 +381 @@
       emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
          &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
-      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar flux
+      !
+      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip
          &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
@@ -379 +390 @@
          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1)
       !
-      CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean
-      CALL iom_put( "qsb_oce", - zqsb )                 ! output downward sensible heat over the ocean
-      CALL iom_put( "qla_oce", - zqla )                 ! output downward latent   heat over the ocean
-      CALL iom_put( "qhc_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
-      CALL iom_put( "qns_oce",   qns  )                 ! output downward non solar heat over the ocean
+#if defined key_lim3
+      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by LIM3)
+      qsr_oce(:,:) = qsr(:,:)
+#endif
+      !
+      IF ( nn_ice == 0 ) THEN
+         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
+         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
+         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
+         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
+         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
+         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
+      ENDIF
       !
       IF(ln_ctl) THEN
@@ -401 +421 @@
  
    
-   SUBROUTINE blk_ice_core(  pst   , pui   , pvi   , palb ,   &
-      &                      p_taui, p_tauj, p_qns , p_qsr,   &
-      &                      p_qla , p_dqns, p_dqla,          &
-      &                      p_tpr , p_spr ,                  &
-      &                      p_fr1 , p_fr2 , cd_grid, pdim  ) 
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE blk_ice_core  ***
+#if defined key_lim2 || defined key_lim3
+   SUBROUTINE blk_ice_core_tau
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_core_tau  ***
       !!
       !! ** Purpose :   provide the surface boundary condition over sea-ice
       !!
-      !! ** Method  :   compute momentum, heat and freshwater exchanged
-      !!                between atmosphere and sea-ice using CORE bulk
-      !!                formulea, ice variables and read atmmospheric fields.
+      !! ** Method  :   compute momentum using CORE bulk
+      !!                formulea, ice variables and read atmospheric fields.
       !!                NB: ice drag coefficient is assumed to be a constant
-      !! 
-      !! caution : the net upward water flux has with mm/day unit
-      !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
-      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pui      ! ice surface velocity (i- and i- components      [m/s]
-      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb     ! ice albedo (all skies)                            [%]
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2]
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s]
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s]
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%]
-      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%]
-      CHARACTER(len=1)          , INTENT(in   ) ::   cd_grid  ! ice grid ( C or B-grid)
-      INTEGER                   , INTENT(in   ) ::   pdim     ! number of ice categories
-      !!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
-      REAL(wp) ::   zst2, zst3
-      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb
-      REAL(wp) ::   zztmp                                        ! temporary variable
-      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f                  ! relative wind module and components at F-point
-      REAL(wp) ::             zwndi_t , zwndj_t                  ! relative wind components at T-point
-      !!
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   z_wnds_t          ! wind speed ( = | U10m - U_ice | ) at T-point
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core')
-      !
-      CALL wrk_alloc( jpi,jpj, z_wnds_t )
-      CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
-
-      ijpl  = pdim                            ! number of ice categories
-
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji, jj    ! dummy loop indices
+      REAL(wp) ::   zcoef_wnorm, zcoef_wnorm2
+      REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point
+      REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_tau')
+      !
       ! local scalars ( place there for vector optimisation purposes)
       zcoef_wnorm  = rhoa * Cice
       zcoef_wnorm2 = rhoa * Cice * 0.5
-      zcoef_dqlw   = 4.0 * 0.95 * Stef
-      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
-      zcoef_dqsb   = rhoa * cpa * Cice
 
 !!gm brutal....
-      z_wnds_t(:,:) = 0.e0
-      p_taui  (:,:) = 0.e0
-      p_tauj  (:,:) = 0.e0
+      utau_ice  (:,:) = 0._wp
+      vtau_ice  (:,:) = 0._wp
+      wndm_ice  (:,:) = 0._wp
 !!gm end
 
-#if defined key_lim3
-      tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1)   ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init
-#endif
       ! ----------------------------------------------------------------------------- !
       !    Wind components and module relative to the moving ocean ( U10m - U_ice )   !
       ! ----------------------------------------------------------------------------- !
-      SELECT CASE( cd_grid )
+      SELECT CASE( cp_ice_msh )
       CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
          !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked)
@@ -484 +460 @@
                ! ... scalar wind at I-point (fld being at T-point)
                zwndi_f = 0.25 * (  sf(jp_wndi)%fnow(ji-1,jj  ,1) + sf(jp_wndi)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pui(ji,jj)
+                  &              + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji  ,jj-1,1)  ) - rn_vfac * u_ice(ji,jj)
                zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   &
-                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * pvi(ji,jj)
+                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj)
                zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )
                ! ... ice stress at I-point
-               p_taui(ji,jj) = zwnorm_f * zwndi_f
-               p_tauj(ji,jj) = zwnorm_f * zwndj_f
+               utau_ice(ji,jj) = zwnorm_f * zwndi_f
+               vtau_ice(ji,jj) = zwnorm_f * zwndj_f
                ! ... scalar wind at T-point (fld being at T-point)
-               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pui(ji,jj+1) + pui(ji+1,jj+1)   &
-                  &                                                    + pui(ji,jj  ) + pui(ji+1,jj  )  )
-               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  pvi(ji,jj+1) + pvi(ji+1,jj+1)   &
-                  &                                                    + pvi(ji,jj  ) + pvi(ji+1,jj  )  )
-               z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   &
+                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  )
+               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   &
+                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  )
+               wndm_ice(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
          END DO
-         CALL lbc_lnk( p_taui  , 'I', -1. )
-         CALL lbc_lnk( p_tauj  , 'I', -1. )
-         CALL lbc_lnk( z_wnds_t, 'T',  1. )
+         CALL lbc_lnk( utau_ice, 'I', -1. )
+         CALL lbc_lnk( vtau_ice, 'I', -1. )
+         CALL lbc_lnk( wndm_ice, 'T',  1. )
          !
       CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean)
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vect. opt.
-               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj  ) + pui(ji,jj) )  )
-               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji  ,jj-1) + pvi(ji,jj) )  )
-               z_wnds_t(ji,jj)  = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
+               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
+               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
             END DO
          END DO
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vect. opt.
-               p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj  ) + z_wnds_t(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) )
-               p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1  ) + z_wnds_t(ji,jj) )                          &
-                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) )
+               utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          &
+                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
+               vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          &
+                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
             END DO
          END DO
-         CALL lbc_lnk( p_taui  , 'U', -1. )
-         CALL lbc_lnk( p_tauj  , 'V', -1. )
-         CALL lbc_lnk( z_wnds_t, 'T',  1. )
+         CALL lbc_lnk( utau_ice, 'U', -1. )
+         CALL lbc_lnk( vtau_ice, 'V', -1. )
+         CALL lbc_lnk( wndm_ice, 'T',  1. )
          !
       END SELECT
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice_core: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
+         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice_core: wndm_ice : ')
+      ENDIF
+
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_tau')
+      
+   END SUBROUTINE blk_ice_core_tau
+
+
+   SUBROUTINE blk_ice_core_flx( ptsu, palb )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE blk_ice_core_flx  ***
+      !!
+      !! ** Purpose :   provide the surface boundary condition over sea-ice
+      !!
+      !! ** Method  :   compute heat and freshwater exchanged
+      !!                between atmosphere and sea-ice using CORE bulk
+      !!                formulea, ice variables and read atmmospheric fields.
+      !! 
+      !! caution : the net upward water flux has with mm/day unit
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu          ! sea ice surface temperature
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb          ! ice albedo (all skies)
+      !!
+      INTEGER  ::   ji, jj, jl    ! dummy loop indices
+      REAL(wp) ::   zst2, zst3
+      REAL(wp) ::   zcoef_dqlw, zcoef_dqla, zcoef_dqsb
+      REAL(wp) ::   zztmp, z1_lsub                               ! temporary variable
+      !!
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw             ! long wave heat flux over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb             ! sensible  heat flux over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw            ! long wave heat sensitivity over ice
+      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb            ! sensible  heat sensitivity over ice
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw       ! evaporation and snw distribution after wind blowing (LIM3)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('blk_ice_core_flx')
+      !
+      CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 
+
+      ! local scalars ( place there for vector optimisation purposes)
+      zcoef_dqlw   = 4.0 * 0.95 * Stef
+      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
+      zcoef_dqsb   = rhoa * cpa * Cice
 
       zztmp = 1. / ( 1. - albo )
       !                                     ! ========================== !
-      DO jl = 1, ijpl                       !  Loop over ice categories  !
+      DO jl = 1, jpl                        !  Loop over ice categories  !
          !                                  ! ========================== !
          DO jj = 1 , jpj
@@ -534 +556 @@
                !      I   Radiative FLUXES   !
                ! ----------------------------!
-               zst2 = pst(ji,jj,jl) * pst(ji,jj,jl)
-               zst3 = pst(ji,jj,jl) * zst2
+               zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
+               zst3 = ptsu(ji,jj,jl) * zst2
                ! Short Wave (sw)
-               p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
+               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
                ! Long  Wave (lw)
-               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * pst(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
+               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
                ! lw sensitivity
                z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3                                               
@@ -549 +571 @@
                ! ... turbulent heat fluxes
                ! Sensible Heat
-               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
+               z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )
                ! Latent Heat
-               p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * z_wnds_t(ji,jj)   &                           
-                  &                         * (  11637800. * EXP( -5897.8 / pst(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
-               ! Latent heat sensitivity for ice (Dqla/Dt)
-               IF( p_qla(ji,jj,jl) > 0._wp ) THEN
-                  p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )
+               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls  * Cice * wndm_ice(ji,jj)   &                           
+                  &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1)  ) )
+              ! Latent heat sensitivity for ice (Dqla/Dt)
+               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
+                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) )
                ELSE
-                  p_dqla(ji,jj,jl) = 0._wp
+                  dqla_ice(ji,jj,jl) = 0._wp
                ENDIF
 
                ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
-               z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj)
+               z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj)
 
                ! ----------------------------!
@@ -567 +589 @@
                ! ----------------------------!
                ! Downward Non Solar flux
-               p_qns (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla (ji,jj,jl)
+               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
                ! Total non solar heat flux sensitivity for ice
-               p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) )
+               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
             END DO
             !
@@ -576 +598 @@
       END DO
       !
+      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s]
+      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
+      CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation
+      CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation
+
+#if defined  key_lim3
+      CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 
+
+      ! --- 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
+
+      ! --- evaporation minus precipitation --- !
+      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
+      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
+      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
+      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
+
+      ! --- heat flux associated with emp --- !
+      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst
+         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
+         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
+         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus )
+
+      ! --- total solar and non solar fluxes --- !
+      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
+      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
+
+      ! --- 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 )
+
+      CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 
+#endif
+
       !--------------------------------------------------------------------
       ! FRACTIONs of net shortwave radiation which is not absorbed in the
@@ -581 +641 @@
       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 )
       !
-      p_fr1(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
-      p_fr2(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
-      !
-      p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s]
-      p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s]
-      CALL iom_put( 'snowpre', p_spr * 86400. )                  ! Snow precipitation
-      CALL iom_put( 'precip' , p_tpr * 86400. )                  ! Total precipitation
+      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
+      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
+      !
       !
       IF(ln_ctl) THEN
-         CALL prt_ctl(tab3d_1=p_qla   , clinfo1=' blk_ice_core: p_qla  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb  : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw  : ', tab3d_2=p_dqla  , clinfo2=' p_dqla : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=p_dqns  , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr   , clinfo2=' p_qsr  : ', kdim=ijpl)
-         CALL prt_ctl(tab3d_1=pst     , clinfo1=' blk_ice_core: pst    : ', tab3d_2=p_qns   , clinfo2=' p_qns  : ', kdim=ijpl)
-         CALL prt_ctl(tab2d_1=p_tpr   , clinfo1=' blk_ice_core: p_tpr  : ', tab2d_2=p_spr   , clinfo2=' p_spr  : ')
-         CALL prt_ctl(tab2d_1=p_taui  , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj  , clinfo2=' p_tauj : ')
-         CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ')
-      ENDIF
-
-      CALL wrk_dealloc( jpi,jpj,   z_wnds_t )
-      CALL wrk_dealloc( jpi,jpj,   pdim, z_qlw, z_qsb, z_dqlw, z_dqsb )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core')
-      !
-   END SUBROUTINE blk_ice_core
+         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice_core: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice_core: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
+         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice_core: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
+         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
+      ENDIF
+
+      CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_core_flx')
+      
+   END SUBROUTINE blk_ice_core_flx
+#endif
 
    SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU,    &

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -21 +21 @@
    USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
+   USE sbcapr
    USE sbcdcy          ! surface boundary condition: diurnal cycle
    USE phycst          ! physical constants
@@ -32 +33 @@
    USE cpl_oasis3      ! OASIS3 coupling
    USE geo2ocean       ! 
-   USE oce   , ONLY : tsn, un, vn
+   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
    USE albedo          !
    USE in_out_manager  ! I/O manager
@@ -40 +41 @@
    USE timing          ! Timing
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE eosbn2
+   USE sbcrnf   , ONLY : l_rnfcpl
 #if defined key_cpl_carbon_cycle
    USE p4zflx, ONLY : oce_co2
@@ -46 +49 @@
    USE ice_domain_size, only: ncat
 #endif
+#if defined key_lim3
+   USE limthd_dh       ! for CALL lim_thd_snwblow
+#endif
+
    IMPLICIT NONE
    PRIVATE
-!EM XIOS-OASIS-MCT compliance
+
    PUBLIC   sbc_cpl_init       ! routine called by sbcmod.F90
    PUBLIC   sbc_cpl_rcv        ! routine called by sbc_ice_lim(_2).F90
@@ -89 +96 @@
    INTEGER, PARAMETER ::   jpr_topm   = 32            ! topmeltn
    INTEGER, PARAMETER ::   jpr_botm   = 33            ! botmeltn
-   INTEGER, PARAMETER ::   jprcv      = 33            ! total number of fields received
-
-   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction 
+   INTEGER, PARAMETER ::   jpr_sflx   = 34            ! salt flux
+   INTEGER, PARAMETER ::   jpr_toce   = 35            ! ocean temperature
+   INTEGER, PARAMETER ::   jpr_soce   = 36            ! ocean salinity
+   INTEGER, PARAMETER ::   jpr_ocx1   = 37            ! ocean current on grid 1
+   INTEGER, PARAMETER ::   jpr_ocy1   = 38            !
+   INTEGER, PARAMETER ::   jpr_ssh    = 39            ! sea surface height
+   INTEGER, PARAMETER ::   jpr_fice   = 40            ! ice fraction          
+   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 ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
    INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
    INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
@@ -106 +122 @@
    INTEGER, PARAMETER ::   jps_ivz1   = 14            !
    INTEGER, PARAMETER ::   jps_co2    = 15
-   INTEGER, PARAMETER ::   jpsnd      = 15            ! total number of fields sended
+   INTEGER, PARAMETER ::   jps_soce   = 16            ! ocean salinity
+   INTEGER, PARAMETER ::   jps_ssh    = 17            ! sea surface height
+   INTEGER, PARAMETER ::   jps_qsroce = 18            ! Qsr above the ocean
+   INTEGER, PARAMETER ::   jps_qnsoce = 19            ! Qns above the ocean
+   INTEGER, PARAMETER ::   jps_oemp   = 20            ! ocean freshwater budget (evap - precip)
+   INTEGER, PARAMETER ::   jps_sflx   = 21            ! salt flux
+   INTEGER, PARAMETER ::   jps_otx1   = 22            ! 2 atmosphere-ocean stress components on grid 1
+   INTEGER, PARAMETER ::   jps_oty1   = 23            ! 
+   INTEGER, PARAMETER ::   jps_rnf    = 24            ! runoffs
+   INTEGER, PARAMETER ::   jps_taum   = 25            ! wind stress module
+   INTEGER, PARAMETER ::   jps_fice2  = 26            ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
+   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
 
    !                                                         !!** namelist namsbc_cpl **
@@ -125 +154 @@
    LOGICAL     ::   ln_usecplmask          !  use a coupling mask file to merge data received from several models
                                            !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
-
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: xcplmask
-
    TYPE ::   DYNARR     
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
@@ -139 +165 @@
 
    !! Substitution
+#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -161 +188 @@
       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,nn_cplmodel) , STAT=ierr(3) )
+      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
       !
       sbc_cpl_alloc = MAXVAL( ierr )
@@ -182 +209 @@
       !!              * initialise the OASIS coupler
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
+      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
       !!
       INTEGER ::   jn   ! dummy loop index
@@ -216 +243 @@
          WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
          WRITE(numout,*)'~~~~~~~~~~~~'
+      ENDIF
+      IF( lwp .AND. ln_cpl ) THEN                        ! control print
          WRITE(numout,*)'  received fields (mutiple ice categogies)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
@@ -359 +388 @@
       srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
+      CASE( 'none'          )       ! nothing to do
       CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
       CASE( 'conservative'  )
@@ -370 +400 @@
       !                                                      !     Runoffs & Calving     !   
       !                                                      ! ------------------------- !
-      srcv(jpr_rnf   )%clname = 'O_Runoff'   ;   IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' )   srcv(jpr_rnf)%laction = .TRUE.
-! This isn't right - really just want ln_rnf_emp changed
-!                                                 IF( TRIM( sn_rcv_rnf%cldes ) == 'climato' )   THEN   ;   ln_rnf = .TRUE.
-!                                                 ELSE                                                 ;   ln_rnf = .FALSE.
-!                                                 ENDIF
+      srcv(jpr_rnf   )%clname = 'O_Runoff'
+      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
+         srcv(jpr_rnf)%laction = .TRUE.
+         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
+         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
+      ENDIF
+      !
       srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
 
@@ -384 +418 @@
       srcv(jpr_qnsmix)%clname = 'O_QnsMix'
       SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
+      CASE( 'none'          )       ! nothing to do
       CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
       CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
@@ -399 +434 @@
       srcv(jpr_qsrmix)%clname = 'O_QsrMix'
       SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
+      CASE( 'none'          )       ! nothing to do
       CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
       CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
@@ -414 +450 @@
       !
       ! non solar sensitivity mandatory for LIM ice model
-      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4) &
+      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4 .AND. nn_components /= jp_iam_sas ) &
          CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
       ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
@@ -447 +483 @@
          srcv(jpr_topm:jpr_botm)%laction = .TRUE.
       ENDIF
-
-      ! Allocate all parts of frcv used for received fields
+      !                                                      ! ------------------------------- !
+      !                                                      !   OPA-SAS coupling - rcv by opa !   
+      !                                                      ! ------------------------------- !
+      srcv(jpr_sflx)%clname = 'O_SFLX'
+      srcv(jpr_fice)%clname = 'RIceFrc'
+      !
+      IF( nn_components == jp_iam_opa ) THEN    ! OPA coupled to SAS via OASIS: force received field by OPA (sent by SAS)
+         srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
+         srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
+         srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
+         srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
+         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
+         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
+         ! Vectors: change of sign at north fold ONLY if on the local grid
+         srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
+         sn_rcv_tau%clvgrd = 'U,V'
+         sn_rcv_tau%clvor = 'local grid'
+         sn_rcv_tau%clvref = 'spherical'
+         sn_rcv_emp%cldes = 'oce only'
+         !
+         IF(lwp) THEN                        ! control print
+            WRITE(numout,*)
+            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
+            WRITE(numout,*)'               OPA component  '
+            WRITE(numout,*)
+            WRITE(numout,*)'  received fields from SAS component '
+            WRITE(numout,*)'                  ice cover '
+            WRITE(numout,*)'                  oce only EMP  '
+            WRITE(numout,*)'                  salt flux  '
+            WRITE(numout,*)'                  mixed oce-ice solar flux  '
+            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
+            WRITE(numout,*)'                  wind stress U,V on local grid and sperical coordinates '
+            WRITE(numout,*)'                  wind stress module'
+            WRITE(numout,*)
+         ENDIF
+      ENDIF
+      !                                                      ! -------------------------------- !
+      !                                                      !   OPA-SAS coupling - rcv by sas  !   
+      !                                                      ! -------------------------------- !
+      srcv(jpr_toce  )%clname = 'I_SSTSST'
+      srcv(jpr_soce  )%clname = 'I_SSSal'
+      srcv(jpr_ocx1  )%clname = 'I_OCurx1'
+      srcv(jpr_ocy1  )%clname = 'I_OCury1'
+      srcv(jpr_ssh   )%clname = 'I_SSHght'
+      srcv(jpr_e3t1st)%clname = 'I_E3T1st'   
+      srcv(jpr_fraqsr)%clname = 'I_FraQsr'   
+      !
+      IF( nn_components == jp_iam_sas ) THEN
+         IF( .NOT. ln_cpl ) srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
+         IF( .NOT. ln_cpl ) srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
+         IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
+         srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
+         srcv( jpr_e3t1st )%laction = lk_vvl
+         srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
+         srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
+         ! Vectors: change of sign at north fold ONLY if on the local grid
+         srcv(jpr_ocx1:jpr_ocy1)%nsgn = -1.
+         ! Change first letter to couple with atmosphere if already coupled OPA
+         ! this is nedeed as each variable name used in the namcouple must be unique:
+         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
+         DO jn = 1, jprcv
+            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
+         END DO
+         !
+         IF(lwp) THEN                        ! control print
+            WRITE(numout,*)
+            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
+            WRITE(numout,*)'               SAS component  '
+            WRITE(numout,*)
+            IF( .NOT. ln_cpl ) THEN
+               WRITE(numout,*)'  received fields from OPA component '
+            ELSE
+               WRITE(numout,*)'  Additional received fields from OPA component : '
+            ENDIF
+            WRITE(numout,*)'               sea surface temperature (Celcius) '
+            WRITE(numout,*)'               sea surface salinity ' 
+            WRITE(numout,*)'               surface currents ' 
+            WRITE(numout,*)'               sea surface height ' 
+            WRITE(numout,*)'               thickness of first ocean T level '        
+            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
+            WRITE(numout,*)
+         ENDIF
+      ENDIF
+      
+      ! =================================================== !
+      ! Allocate all parts of frcv used for received fields !
+      ! =================================================== !
       DO jn = 1, jprcv
          IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
@@ -454 +575 @@
       ! Allocate taum part of frcv which is used even when not received as coupling field
       IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
+      ! Allocate w10m part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
+      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
+      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
+      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
       ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
       IF( k_ice /= 0 ) THEN
@@ -485 +611 @@
       CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
       END SELECT
-     
+           
       !                                                      ! ------------------------- !
       !                                                      !          Albedo           !
@@ -518 +644 @@
          IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
       ENDIF
-
+      
       SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
       CASE( 'none'         )       ! nothing to do
@@ -567 +693 @@
       !                                                      ! ------------------------- !
       ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
+
+      !                                                      ! ------------------------------- !
+      !                                                      !   OPA-SAS coupling - snd by opa !   
+      !                                                      ! ------------------------------- !
+      ssnd(jps_ssh   )%clname = 'O_SSHght' 
+      ssnd(jps_soce  )%clname = 'O_SSSal' 
+      ssnd(jps_e3t1st)%clname = 'O_E3T1st'   
+      ssnd(jps_fraqsr)%clname = 'O_FraQsr'
+      !
+      IF( nn_components == jp_iam_opa ) THEN
+         ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
+         ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
+         ssnd( jps_e3t1st )%laction = lk_vvl
+         ! vector definition: not used but cleaner...
+         ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
+         ssnd(jps_ocy1)%clgrid  = 'V'        !           and           V-point
+         sn_snd_crt%clvgrd = 'U,V'
+         sn_snd_crt%clvor = 'local grid'
+         sn_snd_crt%clvref = 'spherical'
+         !
+         IF(lwp) THEN                        ! control print
+            WRITE(numout,*)
+            WRITE(numout,*)'  sent fields to SAS component '
+            WRITE(numout,*)'               sea surface temperature (T before, Celcius) '
+            WRITE(numout,*)'               sea surface salinity ' 
+            WRITE(numout,*)'               surface currents U,V on local grid and spherical coordinates' 
+            WRITE(numout,*)'               sea surface height ' 
+            WRITE(numout,*)'               thickness of first ocean T level '        
+            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
+            WRITE(numout,*)
+         ENDIF
+      ENDIF
+      !                                                      ! ------------------------------- !
+      !                                                      !   OPA-SAS coupling - snd by sas !   
+      !                                                      ! ------------------------------- !
+      ssnd(jps_sflx  )%clname = 'I_SFLX'     
+      ssnd(jps_fice2 )%clname = 'IIceFrc'
+      ssnd(jps_qsroce)%clname = 'I_QsrOce'   
+      ssnd(jps_qnsoce)%clname = 'I_QnsOce'   
+      ssnd(jps_oemp  )%clname = 'IOEvaMPr' 
+      ssnd(jps_otx1  )%clname = 'I_OTaux1'   
+      ssnd(jps_oty1  )%clname = 'I_OTauy1'   
+      ssnd(jps_rnf   )%clname = 'I_Runoff'   
+      ssnd(jps_taum  )%clname = 'I_TauMod'   
+      !
+      IF( nn_components == jp_iam_sas ) THEN
+         IF( .NOT. ln_cpl ) ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
+         ssnd( (/jps_qsroce, jps_qnsoce, jps_oemp, jps_fice2, jps_sflx, jps_otx1, jps_oty1, jps_taum/) )%laction = .TRUE.
+         !
+         ! Change first letter to couple with atmosphere if already coupled with sea_ice
+         ! this is nedeed as each variable name used in the namcouple must be unique:
+         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
+         DO jn = 1, jpsnd
+            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
+         END DO
+         !
+         IF(lwp) THEN                        ! control print
+            WRITE(numout,*)
+            IF( .NOT. ln_cpl ) THEN
+               WRITE(numout,*)'  sent fields to OPA component '
+            ELSE
+               WRITE(numout,*)'  Additional sent fields to OPA component : '
+            ENDIF
+            WRITE(numout,*)'                  ice cover '
+            WRITE(numout,*)'                  oce only EMP  '
+            WRITE(numout,*)'                  salt flux  '
+            WRITE(numout,*)'                  mixed oce-ice solar flux  '
+            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
+            WRITE(numout,*)'                  wind stress U,V components'
+            WRITE(numout,*)'                  wind stress module'
+         ENDIF
+      ENDIF
+
       !
       ! ================================ !
@@ -572 +771 @@
       ! ================================ !
 
-      CALL cpl_define(jprcv, jpsnd,nn_cplmodel)            
+      CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
+      
       IF (ln_usecplmask) THEN 
          xcplmask(:,:,:) = 0.
@@ -582 +782 @@
          xcplmask(:,:,:) = 1.
       ENDIF
-      !
-      IF( ln_dm2dc .AND. ( cpl_freq( jpr_qsroce ) + cpl_freq( jpr_qsrmix ) /= 86400 ) )   &
+      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
+      !
+      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'S_QsrOce' ) + cpl_freq( 'S_QsrMix' )
+      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
          &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
+      ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
       CALL wrk_dealloc( jpi,jpj, zacs, zaos )
@@ -638 +841 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
-      INTEGER, INTENT(in) ::   k_fsbc   ! frequency of sbc (-> ice model) computation 
-      INTEGER, INTENT(in) ::   k_ice    ! ice management in the sbc (=0/1/2/3)
-      !!
-      LOGICAL ::    llnewtx, llnewtau      ! update wind stress components and module??
+      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
+      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation 
+      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
+
+      !!
+      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
       INTEGER  ::   ji, jj, jn             ! dummy loop indices
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
@@ -650 +854 @@
       REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
       REAL(wp) ::   zzx, zzy               ! temporary variables
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
       !
-      CALL wrk_alloc( jpi,jpj, ztx, zty )
-      !                                                 ! Receive all the atmos. fields (including ice information)
-      isec = ( kt - nit000 ) * NINT( rdttra(1) )             ! date of exchanges
-      DO jn = 1, jprcv                                       ! received fields sent by the atmosphere
-         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask, nrcvinfo(jn) )
+      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
+      !
+      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
+      !
+      !                                                      ! ======================================================= !
+      !                                                      ! Receive all the atmos. fields (including ice information)
+      !                                                      ! ======================================================= !
+      isec = ( kt - nit000 ) * NINT( rdttra(1) )                ! date of exchanges
+      DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
+         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
       END DO
 
@@ -719 +928 @@
          !
       ENDIF
-      
       !                                                      ! ========================= !
       !                                                      !    wind stress module     !   (taum)
@@ -748 +956 @@
          ENDIF
       ENDIF
-      
+      !
       !                                                      ! ========================= !
       !                                                      !      10 m wind speed      !   (wndm)
@@ -761 +969 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi 
-                  wndm(ji,jj) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
+                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
                END DO
             END DO
          ENDIF
-      ELSE
-         IF ( nrcvinfo(jpr_w10m) == OASIS_Rcv ) wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
       ENDIF
 
@@ -773 +979 @@
       IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
          !
-         utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
-         vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
-         taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
+         IF( ln_mixcpl ) THEN
+            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
+            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
+            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
+            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
+         ELSE
+            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
+            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
+            taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
+            wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
+         ENDIF
          CALL iom_put( "taum_oce", taum )   ! output wind stress module
          !  
@@ -781 +995 @@
 
 #if defined key_cpl_carbon_cycle
-      !                                                              ! atmosph. CO2 (ppm)
+      !                                                      ! ================== !
+      !                                                      ! atmosph. CO2 (ppm) !
+      !                                                      ! ================== !
       IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
 #endif
 
+      !  Fields received by SAS when OASIS coupling
+      !  (arrays no more filled at sbcssm stage)
+      !                                                      ! ================== !
+      !                                                      !        SSS         !
+      !                                                      ! ================== !
+      IF( srcv(jpr_soce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
+         sss_m(:,:) = frcv(jpr_soce)%z3(:,:,1)
+         CALL iom_put( 'sss_m', sss_m )
+      ENDIF
+      !                                               
+      !                                                      ! ================== !
+      !                                                      !        SST         !
+      !                                                      ! ================== !
+      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
+         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1)
+         IF( srcv(jpr_soce)%laction .AND. ln_useCT ) THEN    ! make sure that sst_m is the potential temperature
+            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) )
+         ENDIF
+      ENDIF
+      !                                                      ! ================== !
+      !                                                      !        SSH         !
+      !                                                      ! ================== !
+      IF( srcv(jpr_ssh )%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
+         ssh_m(:,:) = frcv(jpr_ssh )%z3(:,:,1)
+         CALL iom_put( 'ssh_m', ssh_m )
+      ENDIF
+      !                                                      ! ================== !
+      !                                                      !  surface currents  !
+      !                                                      ! ================== !
+      IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
+         ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
+         ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
+         CALL iom_put( 'ssu_m', ssu_m )
+      ENDIF
+      IF( srcv(jpr_ocy1)%laction ) THEN
+         ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
+         vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
+         CALL iom_put( 'ssv_m', ssv_m )
+      ENDIF
+      !                                                      ! ======================== !
+      !                                                      !  first T level thickness !
+      !                                                      ! ======================== !
+      IF( srcv(jpr_e3t1st )%laction ) THEN                   ! received by sas in case of opa <-> sas coupling
+         e3t_m(:,:) = frcv(jpr_e3t1st )%z3(:,:,1)
+         CALL iom_put( 'e3t_m', e3t_m(:,:) )
+      ENDIF
+      !                                                      ! ================================ !
+      !                                                      !  fraction of solar net radiation !
+      !                                                      ! ================================ !
+      IF( srcv(jpr_fraqsr)%laction ) THEN                    ! received by sas in case of opa <-> sas coupling
+         frq_m(:,:) = frcv(jpr_fraqsr)%z3(:,:,1)
+         CALL iom_put( 'frq_m', frq_m )
+      ENDIF
+      
       !                                                      ! ========================= !
-      IF( k_ice <= 1 ) THEN                                  !  heat & freshwater fluxes ! (Ocean only case)
+      IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
          !                                                   ! ========================= !
          !
          !                                                       ! total freshwater fluxes over the ocean (emp)
-         SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
-         CASE( 'conservative' )
-            emp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
-         CASE( 'oce only', 'oce and ice' )
-            emp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
-         CASE default
-            CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
-         END SELECT
+         IF( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction ) THEN
+            SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
+            CASE( 'conservative' )
+               zemp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
+            CASE( 'oce only', 'oce and ice' )
+               zemp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
+            CASE default
+               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
+            END SELECT
+         ELSE
+            zemp(:,:) = 0._wp
+         ENDIF
          !
          !                                                        ! runoffs and calving (added in emp)
-         IF( srcv(jpr_rnf)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
-         IF( srcv(jpr_cal)%laction )   emp(:,:) = emp(:,:) - frcv(jpr_cal)%z3(:,:,1)
-         !
-!!gm :  this seems to be internal cooking, not sure to need that in a generic interface 
-!!gm                                       at least should be optional...
-!!         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN     ! add to the total freshwater budget
-!!            ! remove negative runoff
-!!            zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) ) 
-!!            zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
-!!            IF( lk_mpp )   CALL mpp_sum( zcumulpos )   ! sum over the global domain
-!!            IF( lk_mpp )   CALL mpp_sum( zcumulneg ) 
-!!            IF( zcumulpos /= 0. ) THEN                 ! distribute negative runoff on positive runoff grid points
-!!               zcumulneg = 1.e0 + zcumulneg / zcumulpos
-!!               frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
-!!            ENDIF     
-!!            ! add runoff to e-p 
-!!            emp(:,:) = emp(:,:) - frcv(jpr_rnf)%z3(:,:,1)
-!!         ENDIF
-!!gm  end of internal cooking
+         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
+         
+         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
+         ELSE                   ;   emp(:,:) =                              zemp(:,:)
+         ENDIF
          !
          !                                                       ! non solar heat flux over the ocean (qns)
-         IF( srcv(jpr_qnsoce)%laction )   qns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
-         IF( srcv(jpr_qnsmix)%laction )   qns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
+         IF(      srcv(jpr_qnsoce)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
+         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
+         ELSE                                       ;   zqns(:,:) = 0._wp
+         END IF
          ! update qns over the free ocean with:
-         qns(:,:) =  qns(:,:) - emp(:,:) * sst_m(:,:) * rcp            ! remove heat content due to mass flux (assumed to be at SST)
-         IF( srcv(jpr_snow  )%laction )   THEN
-              qns(:,:) = qns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
+         IF( nn_components /= jp_iam_opa ) THEN
+            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
+            IF( srcv(jpr_snow  )%laction ) THEN
+               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
+            ENDIF
+         ENDIF
+         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
+         ELSE                   ;   qns(:,:) =                              zqns(:,:)
          ENDIF
 
          !                                                       ! solar flux over the ocean          (qsr)
-         IF( srcv(jpr_qsroce)%laction )   qsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
-         IF( srcv(jpr_qsrmix)%laction )   qsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
-         IF( ln_dm2dc )   qsr(:,:) = sbc_dcy( qsr )                           ! modify qsr to include the diurnal cycle
+         IF     ( srcv(jpr_qsroce)%laction ) THEN   ;   zqsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
+         ELSE IF( srcv(jpr_qsrmix)%laction ) then   ;   zqsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
+         ELSE                                       ;   zqsr(:,:) = 0._wp
+         ENDIF
+         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
+         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
+         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
+         ENDIF
          !
-  
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj, ztx, zty )
+         ! salt flux over the ocean (received by opa in case of opa <-> sas coupling)
+         IF( srcv(jpr_sflx )%laction )   sfx(:,:) = frcv(jpr_sflx  )%z3(:,:,1)
+         ! Ice cover  (received by opa in case of opa <-> sas coupling)
+         IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
+         !
+
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
@@ -934 +1212 @@
             !
          ENDIF
-
          !                                                      ! ======================= !
          !                                                      !     put on ice grid     !
@@ -1056 +1333 @@
    
 
-   SUBROUTINE sbc_cpl_ice_flx( p_frld  , palbi   , psst    , pist    )
+   SUBROUTINE sbc_cpl_ice_flx( p_frld, palbi, psst, pist )
       !!----------------------------------------------------------------------
       !!             ***  ROUTINE sbc_cpl_ice_flx  ***
@@ -1098 +1375 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
       ! optional arguments, used only in 'mixed oce-ice' case
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi   ! all skies ice albedo 
-      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst    ! sea surface temperature     [Celsius]
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist    ! ice surface temperature     [Kelvin]
-      !
-      INTEGER ::   jl   ! dummy loop index
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zcptn, ztmp, zicefr
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo 
+      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
+      !
+      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, zqprec_ice, zqemp_oce ! for LIM3
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
       !
-      CALL wrk_alloc( jpi,jpj, zcptn, ztmp, zicefr )
-
+      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 )
+
+      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
       zicefr(:,:) = 1.- p_frld(:,:)
       zcptn(:,:) = rcp * sst_m(:,:)
@@ -1117 +1399 @@
       !                                                      ! ========================= !
       !
-      !                                                           ! total Precipitations - total Evaporation (emp_tot)
-      !                                                           ! solid precipitation  - sublimation       (emp_ice)
-      !                                                           ! solid Precipitation                      (sprecip)
+      !                                                           ! total Precipitation - total Evaporation (emp_tot)
+      !                                                           ! solid precipitation - sublimation       (emp_ice)
+      !                                                           ! solid Precipitation                     (sprecip)
+      !                                                           ! liquid + solid Precipitation            (tprecip)
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
       CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
-         sprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                 ! May need to ensure positive here
-         tprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + sprecip (:,:) ! May need to ensure positive here
-         emp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - tprecip(:,:)
-         emp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
+         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 
          IF( iom_use('hflx_rain_cea') )   &
@@ -1136 +1419 @@
             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
-         emp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
-         emp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
-         sprecip(:,:) = - frcv(jpr_semp)%z3(:,:,1) + frcv(jpr_ievp)%z3(:,:,1)
+         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
+         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
+         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( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
+      IF( srcv(jpr_cal)%laction ) THEN 
+         zemp_tot(:,:) = zemp_tot(:,:) - 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(:,:)
+         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
+         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
+      ELSE
+         emp_tot(:,:) =                                  zemp_tot(:,:)
+         emp_ice(:,:) =                                  zemp_ice(:,:)
+         sprecip(:,:) =                                  zsprecip(:,:)
+         tprecip(:,:) =                                  ztprecip(:,:)
+      ENDIF
 
          CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
@@ -1146 +1452 @@
       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)
-      !   
-      !                                                           ! runoffs and calving (put in emp_tot)
-      IF( srcv(jpr_rnf)%laction ) THEN 
-         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)
-            CALL iom_put( 'runoffs'      , frcv(jpr_rnf)%z3(:,:,1)              )   ! rivers
-         IF( iom_use('hflx_rnf_cea') )   &
-            CALL iom_put( 'hflx_rnf_cea' , frcv(jpr_rnf)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from rivers
-      ENDIF
-      IF( srcv(jpr_cal)%laction ) THEN 
-         emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
-         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
-      ENDIF
-      !
-!!gm :  this seems to be internal cooking, not sure to need that in a generic interface 
-!!gm                                       at least should be optional...
-!!       ! remove negative runoff                            ! sum over the global domain
-!!       zcumulpos = SUM( MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) ) 
-!!       zcumulneg = SUM( MIN( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * e1t(:,:) * e2t(:,:) * tmask_i(:,:) )
-!!       IF( lk_mpp )   CALL mpp_sum( zcumulpos )
-!!       IF( lk_mpp )   CALL mpp_sum( zcumulneg ) 
-!!       IF( zcumulpos /= 0. ) THEN                          ! distribute negative runoff on positive runoff grid points
-!!          zcumulneg = 1.e0 + zcumulneg / zcumulpos
-!!          frcv(jpr_rnf)%z3(:,:,1) = MAX( frcv(jpr_rnf)%z3(:,:,1), 0.e0 ) * zcumulneg
-!!       ENDIF     
-!!       emp_tot(:,:) = emp_tot(:,:) - frcv(jpr_rnf)%z3(:,:,1)   ! add runoff to e-p 
-!!
-!!gm  end of internal cooking
 
       !                                                      ! ========================= !
@@ -1180 +1457 @@
       !                                                      ! ========================= !
       CASE( 'oce only' )                                     ! the required field is directly provided
-         qns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
+         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
       CASE( 'conservative' )                                      ! the required fields are directly provided
-         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
+         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
-            qns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
+            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
          ELSE
             ! Set all category values equal for the moment
             DO jl=1,jpl
-               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
             ENDDO
          ENDIF
       CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
-         qns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
+         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
             DO jl=1,jpl
-               qns_tot(:,:   ) = qns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
-               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
+               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
             ENDDO
          ELSE
             qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl=1,jpl
-               qns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
             ENDDO
          ENDIF
       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 **
-         qns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
-         qns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
+         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(:,:) ) )
       END SELECT
-      ztmp(:,:) = p_frld(:,:) * sprecip(:,:) * lfus
-      qns_tot(:,:) = qns_tot(:,:)                         &            ! qns_tot update over free ocean with:
-         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
-         &          - (  emp_tot(:,:)                     &            ! remove the heat content of mass flux (assumed to be at SST)
-         &             - emp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
-      IF( iom_use('hflx_snow_cea') )   &
-         CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
 !!gm
-!!    currently it is taken into account in leads budget but not in the qns_tot, and thus not in 
+!!    currently it is taken into account in leads budget but not in the zqns_tot, and thus not in 
 !!    the flux that enter the ocean....
 !!    moreover 1 - it is not diagnose anywhere.... 
@@ -1228 +1499 @@
       IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting 
          ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting 
-         qns_tot(:,:) = qns_tot(:,:) - ztmp(:,:)
+         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(:,:)
+
+      ! --- non solar flux over ocean --- !
+      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
+      zqns_oce = 0._wp
+      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
+
+      ! --- heat flux associated with emp --- !
+      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) --- !
+      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
+
+      ! --- total non solar flux --- !
+      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
+
+      ! --- in case both coupled/forced are active, we must mix values --- ! 
+      IF( ln_mixcpl ) THEN
+         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
+         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
+         DO jl=1,jpl
+            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
+         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)
+      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 ) 
+
+#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
+         &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
+         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
+
+     IF( ln_mixcpl ) THEN
+         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
+         DO jl=1,jpl
+            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
+         ENDDO
+      ELSE
+         qns_tot(:,:  ) = zqns_tot(:,:  )
+         qns_ice(:,:,:) = zqns_ice(:,:,:)
+      ENDIF
+
+#endif
 
       !                                                      ! ========================= !
@@ -1237 +1586 @@
       !                                                      ! ========================= !
       CASE( 'oce only' )
-         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
+         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
       CASE( 'conservative' )
-         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
+         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
-            qsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
+            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
          ELSE
             ! Set all category values equal for the moment
             DO jl=1,jpl
-               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             ENDDO
          ENDIF
-         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
-         qsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
+         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
+         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
       CASE( 'oce and ice' )
-         qsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
+         zqsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
             DO jl=1,jpl
-               qsr_tot(:,:   ) = qsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
-               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
+               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
             ENDDO
          ELSE
             qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
             DO jl=1,jpl
-               qsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
+               zqsr_tot(:,:   ) = zqsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
+               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             ENDDO
          ENDIF
       CASE( 'mixed oce-ice' )
-         qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
+         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
 ! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
 !       Create solar heat flux over ice using incoming solar heat flux and albedos
 !       ( see OASIS3 user guide, 5th edition, p39 )
-         qsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
+         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
             &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
             &                     + palbi         (:,:,1) * zicefr(:,:) ) )
       END SELECT
-      IF( ln_dm2dc ) THEN   ! modify qsr to include the diurnal cycle
-         qsr_tot(:,:  ) = sbc_dcy( qsr_tot(:,:  ) )
+      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
+         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
          DO jl=1,jpl
-            qsr_ice(:,:,jl) = sbc_dcy( qsr_ice(:,:,jl) )
+            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
          ENDDO
+      ENDIF
+
+      IF( ln_mixcpl ) THEN
+         qsr_tot(:,:) = qsr(:,:) * p_frld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
+         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
+         DO jl=1,jpl
+            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
+         ENDDO
+      ELSE
+         qsr_tot(:,:  ) = zqsr_tot(:,:  )
+         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
       ENDIF
 
@@ -1284 +1645 @@
       CASE ('coupled')
          IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
-            dqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
+            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
          ELSE
             ! Set all category values equal for the moment
             DO jl=1,jpl
-               dqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
+               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
             ENDDO
          ENDIF
       END SELECT
-
+      
+      IF( ln_mixcpl ) THEN
+         DO jl=1,jpl
+            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
+         ENDDO
+      ELSE
+         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
+      ENDIF
+      
       !                                                      ! ========================= !
       SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
@@ -1308 +1677 @@
       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
 
-      CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
+      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 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
@@ -1328 +1698 @@
       INTEGER ::   ji, jj, jl   ! dummy loop indices
       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(:,:,:) ::   ztmp3, ztmp4   
@@ -1344 +1715 @@
       !                                                      ! ------------------------- !
       IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
-         SELECT CASE( sn_snd_temp%cldes)
-         CASE( 'oce only'             )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0
-         CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:)   
-            SELECT CASE( sn_snd_temp%clcat )
-            CASE( 'yes' )   
-               ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
-            CASE( 'no' )
-               ztmp3(:,:,:) = 0.0
+         
+         IF ( nn_components == jp_iam_opa ) THEN
+            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of ln_useCT on the received part
+         ELSE
+            ! we must send the surface potential temperature 
+            IF( ln_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
+            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
+            ENDIF
+            !
+            SELECT CASE( sn_snd_temp%cldes)
+            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
+            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
+               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(:,:) 
                DO jl=1,jpl
-                  ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
+                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
                ENDDO
-            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
+            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
             END SELECT
-         CASE( 'mixed oce-ice'        )   
-            ztmp1(:,:) = ( tsn(:,:,1,jp_tem) + rt0 ) * zfr_l(:,:) 
-            DO jl=1,jpl
-               ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
-            ENDDO
-         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
-         END SELECT
+         ENDIF
          IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
          IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
@@ -1385 +1766 @@
       !                                                      !  Ice fraction & Thickness ! 
       !                                                      ! ------------------------- !
-      ! Send ice fraction field 
+      ! Send ice fraction field to atmosphere
       IF( ssnd(jps_fice)%laction ) THEN
          SELECT CASE( sn_snd_thick%clcat )
@@ -1392 +1773 @@
          CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
          END SELECT
-         CALL cpl_snd( jps_fice, isec, ztmp3, info )
+         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
+      ENDIF
+      
+      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
+      IF( ssnd(jps_fice2)%laction ) THEN
+         ztmp3(:,:,1) = fr_i(:,:)
+         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
       ENDIF
 
@@ -1440 +1827 @@
          !                                                              i-1  i   i
          !                                                               i      i+1 (for I)
-         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) ) 
-               END DO
-            END DO
-         CASE( 'weighted oce and ice' )   
-            SELECT CASE ( cp_ice_msh )
-            CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
+         IF( nn_components == jp_iam_opa ) THEN
+            zotx1(:,:) = un(:,:,1)  
+            zoty1(:,:) = vn(:,:,1)  
+         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) ) * zfr_l(ji,jj)  
-                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
-                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
-                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+                     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
-            CASE( 'I' )                      ! Ocean on C grid, Ice on I-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)  
-                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)  
-                     zitx1(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)
-                     zity1(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)
+            CASE( 'weighted oce and ice' )   
+               SELECT CASE ( cp_ice_msh )
+               CASE( 'C' )                      ! Ocean and Ice on 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) ) * zfr_l(ji,jj)  
+                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
+                        zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
+                        zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+                     END DO
                   END DO
-               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)  
-                     zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)  
-                     zitx1(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)
-                     zity1(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)
+               CASE( 'I' )                      ! Ocean on C grid, Ice on I-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)  
+                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)  
+                        zitx1(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)
+                        zity1(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
+               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)  
+                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)  
+                        zitx1(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)
+                        zity1(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 SELECT
+               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
+            CASE( 'mixed oce-ice'        )
+               SELECT CASE ( cp_ice_msh )
+               CASE( 'C' )                      ! Ocean and Ice on 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) ) * zfr_l(ji,jj)   &
+                           &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
+                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
+                           &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
+                     END DO
+                  END DO
+               CASE( 'I' )                      ! Ocean on C grid, Ice on I-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)
+                     END DO
+                  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)
+                     END DO
+                  END DO
+               END SELECT
             END SELECT
-            CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
-         CASE( 'mixed oce-ice'        )
-            SELECT CASE ( cp_ice_msh )
-            CASE( 'C' )                      ! Ocean and Ice on 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) ) * zfr_l(ji,jj)   &
-                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
-                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
-                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
-                  END DO
-               END DO
-            CASE( 'I' )                      ! Ocean on C grid, Ice on I-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)
-                  END DO
-               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)
-                  END DO
-               END DO
-            END SELECT
-         END SELECT
-         CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
+            CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
+            !
+         ENDIF
          !
          !
@@ -1558 +1951 @@
       ENDIF
       !
+      !
+      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
+      !                                                        ! SSH
+      IF( ssnd(jps_ssh )%laction )  THEN
+         !                          ! removed inverse barometer ssh when Patm
+         !                          forcing is used (for sea-ice dynamics)
+         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
+         ENDIF
+         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
+
+      ENDIF
+      !                                                        ! SSS
+      IF( ssnd(jps_soce  )%laction )  THEN
+         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
+      ENDIF
+      !                                                        ! first T level thickness 
+      IF( ssnd(jps_e3t1st )%laction )  THEN
+         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( fse3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
+      ENDIF
+      !                                                        ! Qsr fraction
+      IF( ssnd(jps_fraqsr)%laction )  THEN
+         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
+      ENDIF
+      !
+      !  Fields sent by SAS to OPA when OASIS coupling
+      !                                                        ! Solar heat flux
+      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
+      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
+      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
+      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
+      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
+      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
+      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 )
+
       CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
       CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

@@ -138 +138 @@
          IF      ( ksbc == jp_flx ) THEN
             CALL cice_sbc_force(kt)
-         ELSE IF ( ksbc == jp_cpl ) THEN
+         ELSE IF ( ksbc == jp_purecpl ) THEN
             CALL sbc_cpl_ice_flx( 1.0-fr_i  )
          ENDIF
@@ -146 +146 @@
          CALL cice_sbc_out ( kt, ksbc )
 
-         IF ( ksbc == jp_cpl )  CALL cice_sbc_hadgam(kt+1)
+         IF ( ksbc == jp_purecpl )  CALL cice_sbc_hadgam(kt+1)
 
       ENDIF                                          ! End sea-ice time step only
@@ -187 +187 @@
 
 ! Do some CICE consistency checks
-      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_cpl) ) THEN
+      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          IF ( calc_strair .OR. calc_Tsfc ) THEN
             CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
@@ -212 +212 @@
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
-      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_cpl) ) THEN
+      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
@@ -319 +319 @@
 ! forced and coupled case 
 
-      IF ( (ksbc == jp_flx).OR.(ksbc == jp_cpl) ) THEN
+      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
 
          ztmpn(:,:,:)=0.0
@@ -587 +587 @@
       ELSE IF (ksbc == jp_core) THEN
          emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)        
-      ELSE IF (ksbc == jp_cpl) THEN
+      ELSE IF (ksbc == jp_purecpl) THEN
 ! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above) 
 ! This is currently as required with the coupling fields from the UM atmosphere
@@ -623 +623 @@
       ENDIF
 ! Take into account snow melting except for fully coupled when already in qns_tot
-      IF (ksbc == jp_cpl) THEN
+      IF (ksbc == jp_purecpl) THEN
          qsr(:,:)= qsr_tot(:,:)
          qns(:,:)= qns_tot(:,:)
@@ -658 +658 @@
 
       CALL cice2nemo(aice,fr_i,'T', 1. )
-      IF ( (ksbc == jp_flx).OR.(ksbc == jp_cpl) ) THEN
+      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
          DO jl=1,ncat
             CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90

@@ -105 +105 @@
          fr_i(:,:) = eos_fzp( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius]
 
-         IF( lk_cpl )   a_i(:,:,1) = fr_i(:,:)         
+         IF( ln_cpl )   a_i(:,:,1) = fr_i(:,:)         
 
          ! Flux and ice fraction computation

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

@@ -110 +110 @@
       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 
       !!----------------------------------------------------------------------
 
@@ -115 +116 @@
 
       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only
+
          !-----------------------!                                           
          ! --- Bulk Formulae --- !                                           
@@ -124 +126 @@
          t_bo(:,:) = ( eos_fzp( sss_m ) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )  
          !                                                                                      
-         ! Ice albedo
-         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
-         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
-
-         ! CORE and COUPLED bulk formulations
-         SELECT CASE( kblk )
-         CASE( jp_core , jp_cpl )
-
-            ! albedo depends on cloud fraction because of non-linear spectral effects
-            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
-            ! (zalb_ice) is computed within the bulk routine
-            
-         END SELECT
+!!clem         ! Ice albedo
+!!clem         CALL wrk_@lloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+!!clem         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
+!!
+!!         ! CORE and COUPLED bulk formulations
+!!         SELECT CASE( kblk )
+!!         CASE( jp_core , jp_purecpl )
+!!            ! albedo depends on cloud fraction because of non-linear spectral effects
+!!            zalb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+!!            ! In CLIO the cloud fraction is read in the climatology and the all-sky albedo 
+!!            ! (zalb_ice) is computed within the bulk routine
+!!clem         END SELECT
          
          ! Mask sea ice surface temperature (set to rt0 over land)
@@ -154 +154 @@
          SELECT CASE( kblk )
          CASE( jp_clio )                                       ! CLIO bulk formulation
-            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               &
-               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
-               &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
-               &                      tprecip    , sprecip    ,                           &
-               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
-            !         
-            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+!!clem            CALL blk_ice_clio( t_su , zalb_cs    , zalb_os    , zalb_ice  ,               &
+!!               &                      utau_ice   , vtau_ice   , qns_ice   , qsr_ice   ,   &
+!!               &                      qla_ice    , dqns_ice   , dqla_ice  ,               &
+!!               &                      tprecip    , sprecip    ,                           &
+!!               &                      fr1_i0     , fr2_i0     , cp_ice_msh, jpl  )
+!!            !         
+!!            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
+!!               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+            CALL blk_ice_clio_tau
 
          CASE( jp_core )                                       ! CORE bulk formulation
-            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               &
-               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
-               &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
-               &                      tprecip   , sprecip   ,                            &
-               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
-               !
-            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+!!clem            CALL blk_ice_core( t_su , u_ice     , v_ice     , zalb_ice   ,               &
+!!clem               &                      utau_ice  , vtau_ice  , qns_ice    , qsr_ice   ,   &
+!!clem               &                      qla_ice   , dqns_ice  , dqla_ice   ,               &
+!!clem               &                      tprecip   , sprecip   ,                            &
+!!clem               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
+!!clem            IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
+!!clem               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+            CALL blk_ice_core_tau
             !
-         CASE ( jp_cpl )
+         CASE ( jp_purecpl )
             
             CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
@@ -179 +180 @@
          END SELECT
          
-         !------------------------------!
-         ! --- LIM-3 main time-step --- !
-         !------------------------------!
+         IF( ln_mixcpl) THEN
+            CALL wrk_alloc( jpi,jpj    , zutau_ice, zvtau_ice)
+            CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
+            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+            CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice)
+         ENDIF
+
+         !                                           !----------------------!
+         !                                           ! LIM-3  time-stepping !
+         !                                           !----------------------!
+         ! 
          numit = numit + nn_fsbc                     ! Ice model time step
          !                                                   
@@ -220 +230 @@
          phicif(:,:)  = vt_i(:,:)
          
+         ! Ice albedo
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
+ 
          SELECT CASE( kblk )
-         CASE ( jp_cpl )
-            CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=zalb_ice, psst=sst_m, pist=t_su    )
+         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, qns_ice   , qsr_ice   ,    &
+!              &                      qla_ice, dqns_ice   , dqla_ice  , tprecip, sprecip    ,  &
+!              &                      fr1_i0     , fr2_i0     , jpl  )
+!           !         
+            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 )
+
+         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 )
+
+         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 )
             IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-               &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
+               &                                           dqns_ice, evap_ice, devap_ice, nn_limflx )
             ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux)
-            qla_ice  (:,:,:) = 0._wp
-            dqla_ice (:,:,:) = 0._wp
+            evap_ice  (:,:,:) = 0._wp
+            devap_ice (:,:,:) = 0._wp
+
          END SELECT
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
+
          !
          CALL lim_thd( kt )                         ! Ice thermodynamics 
@@ -247 +286 @@
          IF( ln_icectl )  CALL lim_ctl( kt )        ! alerts in case of model crash
          !
-         CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs, zalb_ice )
          !
       ENDIF   ! End sea-ice time step only
@@ -476 +514 @@
    
    SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice,   &
-         &                          pdqn_ice, pqla_ice, pdql_ice, k_limflx )
+         &                          pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE ice_lim_flx  ***
@@ -494 +532 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqsr_ice   ! net solar flux
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdqn_ice   ! non solar flux sensitivity
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pqla_ice   ! latent heat flux
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdql_ice   ! latent heat flux sensitivity
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pevap_ice  ! sublimation
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdevap_ice ! sublimation sensitivity
       !
       INTEGER  ::   jl      ! dummy loop index
@@ -504 +542 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m   ! Mean solar heat flux over all categories
       REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m   ! Mean non solar heat flux over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_m   ! Mean latent heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m  ! Mean sublimation over all categories
       REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m   ! Mean d(qns)/dT over all categories
-      REAL(wp), POINTER, DIMENSION(:,:) :: z_dql_m   ! Mean d(qla)/dT over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
       !!----------------------------------------------------------------------
 
@@ -514 +552 @@
       SELECT CASE( k_limflx )                              !==  averaged on all ice categories  ==!
       CASE( 0 , 1 )
-         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
+         CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
          !
          z_qns_m(:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
          z_qsr_m(:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
          z_dqn_m(:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
-         z_qla_m(:,:) = fice_ice_ave ( pqla_ice (:,:,:) )
-         z_dql_m(:,:) = fice_ice_ave ( pdql_ice (:,:,:) )
+         z_evap_m(:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
+         z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
          DO jl = 1, jpl
             pdqn_ice(:,:,jl) = z_dqn_m(:,:)
-            pdql_ice(:,:,jl) = z_dql_m(:,:)
+            pdevap_ice(:,:,jl) = z_devap_m(:,:)
          END DO
          !
@@ -529 +567 @@
             pqns_ice(:,:,jl) = z_qns_m(:,:)
             pqsr_ice(:,:,jl) = z_qsr_m(:,:)
-            pqla_ice(:,:,jl) = z_qla_m(:,:)
+            pevap_ice(:,:,jl) = z_evap_m(:,:)
          END DO
          !
-         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_qla_m, z_dqn_m, z_dql_m)
+         CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
       END SELECT
 
@@ -543 +581 @@
          DO jl = 1, jpl
             pqns_ice(:,:,jl) = pqns_ice(:,:,jl) + pdqn_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
-            pqla_ice(:,:,jl) = pqla_ice(:,:,jl) + pdql_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
+            pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * (ptn_ice(:,:,jl) - ztem_m(:,:))
             pqsr_ice(:,:,jl) = pqsr_ice(:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) ) 
          END DO
@@ -593 +631 @@
       wfx_spr(:,:) = 0._wp   ;   
       
-      hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
       hfx_thd(:,:) = 0._wp   ;   
       hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
@@ -610 +647 @@
       
    END SUBROUTINE sbc_lim_diag0
-      
+
+     
    FUNCTION fice_cell_ave ( ptab )
       !!--------------------------------------------------------------------------

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90

@@ -101 +101 @@
       REAL(wp), DIMENSION(:,:,:), POINTER :: zalb_ice  ! mean ice albedo
       REAL(wp), DIMENSION(:,:,:), POINTER :: zsist     ! ice surface temperature (K)
+      REAL(wp), DIMENSION(:,:  ), POINTER :: zutau_ice, zvtau_ice 
       !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi,jpj,1, zalb_os, zalb_cs, zalb_ice, zsist )
 
       IF( kt == nit000 ) THEN
@@ -124 +123 @@
          &*Agrif_PArent(nn_fsbc)/REAL(nn_fsbc)) + 1
 # endif
+
+         CALL wrk_alloc( jpi,jpj  , zutau_ice, zvtau_ice)
+         CALL wrk_alloc( jpi,jpj,1, zalb_os, zalb_cs, zalb_ice, zsist )
+
          !  Bulk Formulea !
          !----------------!
@@ -132 +135 @@
                DO ji = 2, jpi   ! NO vector opt. possible
                   u_oce(ji,jj) = 0.5_wp * ( ssu_m(ji-1,jj  ) * umask(ji-1,jj  ,1) &
-                     &           + ssu_m(ji-1,jj-1) * umask(ji-1,jj-1,1) ) * tmu(ji,jj)
+                     &                    + ssu_m(ji-1,jj-1) * umask(ji-1,jj-1,1) ) * tmu(ji,jj)
                   v_oce(ji,jj) = 0.5_wp * ( ssv_m(ji  ,jj-1) * vmask(ji  ,jj-1,1) &
-                     &           + ssv_m(ji-1,jj-1) * vmask(ji-1,jj-1,1) ) * tmu(ji,jj)
+                     &                    + ssv_m(ji-1,jj-1) * vmask(ji-1,jj-1,1) ) * tmu(ji,jj)
                END DO
             END DO
@@ -158 +161 @@
 
          SELECT CASE( ksbc )
-         CASE( jp_core , jp_cpl )   ! CORE and COUPLED bulk formulations
+         CASE( jp_core , jp_purecpl )   ! CORE and COUPLED bulk formulations
 
             ! albedo depends on cloud fraction because of non-linear spectral effects
@@ -182 +185 @@
          SELECT CASE( ksbc )
          CASE( jp_clio )           ! CLIO bulk formulation
-            CALL blk_ice_clio( zsist, zalb_cs    , zalb_os    , zalb_ice   ,            &
-               &                      utau_ice   , vtau_ice   , qns_ice    , qsr_ice,   &
-               &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
-               &                      tprecip    , sprecip    ,                         &
-               &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
+!           CALL blk_ice_clio( zsist, zalb_cs    , zalb_os    , zalb_ice   ,            &
+!              &                      utau_ice   , vtau_ice   , qns_ice    , qsr_ice,   &
+!              &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
+!              &                      tprecip    , sprecip    ,                         &
+!              &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
+            CALL blk_ice_clio_tau
+            CALL blk_ice_clio_flx( zsist, zalb_cs, zalb_os, zalb_ice )
 
          CASE( jp_core )           ! CORE bulk formulation
-            CALL blk_ice_core( zsist, u_ice      , v_ice      , zalb_ice   ,            &
-               &                      utau_ice   , vtau_ice   , qns_ice    , qsr_ice,   &
-               &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
-               &                      tprecip    , sprecip    ,                         &
-               &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
-
-         CASE( jp_cpl )            ! Coupled formulation : atmosphere-ice stress only (fluxes provided after ice dynamics)
+            CALL blk_ice_core_tau
+            CALL blk_ice_core_flx( zsist, zalb_ice )
+
+         CASE( jp_purecpl )            ! Coupled formulation : atmosphere-ice stress only (fluxes provided after ice dynamics)
             CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
          END SELECT
+         
+         IF( ln_mixcpl) THEN
+            CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
+            utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+            vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
+         ENDIF
 
          CALL iom_put( 'utau_ice', utau_ice )     ! Wind stress over ice along i-axis at I-point
@@ -227 +235 @@
          END IF
          !                                             ! Ice surface fluxes in coupled mode 
-         IF( ksbc == jp_cpl )   THEN
+         IF( ln_cpl ) THEN   ! pure coupled and mixed forced-coupled configurations
             a_i(:,:,1)=fr_i
             CALL sbc_cpl_ice_flx( frld,                                              &
@@ -249 +257 @@
          IF( .NOT. Agrif_Root() )   CALL agrif_update_lim2( kt )
 # endif
+         !
+         CALL wrk_dealloc( jpi,jpj  , zutau_ice, zvtau_ice)
+         CALL wrk_dealloc( jpi,jpj,1, zalb_os, zalb_cs, zalb_ice, zsist )
          !
       ENDIF                                    ! End sea-ice time step only
@@ -260 +271 @@
       IF( ln_limdyn    )   CALL lim_sbc_tau_2( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
       !
-      CALL wrk_dealloc( jpi,jpj,1, zalb_os, zalb_cs, zalb_ice, zsist )
-      !
    END SUBROUTINE sbc_ice_lim_2
 

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

@@ -39 +39 @@
    USE sbcice_cice      ! surface boundary condition: CICE    sea-ice model
    USE sbccpl           ! surface boundary condition: coupled florulation
+   USE cpl_oasis3       ! OASIS routines for coupling
    USE sbcssr           ! surface boundary condition: sea surface restoring
    USE sbcrnf           ! surface boundary condition: runoffs
@@ -84 +85 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core,           &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
-         &             ln_ssr    ,  nn_isf , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, nn_limflx
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
+         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
+         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
       INTEGER  ::   ios
+      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
+      LOGICAL  ::   ll_purecpl
       !!----------------------------------------------------------------------
 
@@ -115 +119 @@
           nn_ice      =   0
       ENDIF
-     
+
       IF(lwp) THEN               ! Control print
          WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
@@ -125 +129 @@
          WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
-         WRITE(numout,*) '              coupled    formulation (T if key_oasis3)   lk_cpl      = ', lk_cpl
+         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
+         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
+         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
+         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
          WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
          WRITE(numout,*) '           Misc. options of sbc : '
@@ -152 +159 @@
       END SELECT
       !
+      IF ( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
+      IF ( nn_components == jp_iam_opa .AND. ln_cpl )   &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
+      IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      IF ( ln_cpl .AND. .NOT. lk_oasis )    &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
+      IF( ln_mixcpl .AND. .NOT. lk_oasis )    &
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' )
+      IF( ln_mixcpl .AND. .NOT. ln_cpl )    &
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' )
+      IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    &
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' )
+
       !                              ! allocate sbc arrays
       IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
@@ -170 +192 @@
          fwfisf_b(:,:) = 0.0_wp
       END IF
-      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero
+      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero
 
       sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
@@ -180 +202 @@
 
       !                                            ! restartability   
-      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
-          MOD( nstock             , nn_fsbc) /= 0 ) THEN 
-         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
-            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
-         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
-      ENDIF
-      !
-      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
-         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
-      !
-      IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. lk_cpl ) )   &
+      IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. ln_cpl ) )   &
          &   CALL ctl_stop( 'LIM sea-ice model requires a bulk formulation or coupled configuration' )
-      IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. lk_cpl ) )   &
-         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or lk_cpl' )
+      IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. ln_cpl ) )   &
+         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or ln_cpl' )
       IF( nn_ice == 4 .AND. lk_agrif )   &
          &   CALL ctl_stop( 'CICE sea-ice model not currently available with AGRIF' )
@@ -200 +212 @@
       IF( ( nn_ice /= 3 ) .AND. ( nn_limflx >= 0 ) )   &
          &   WRITE(numout,*) 'The nn_limflx>=0 option has no effect if sea ice model is not LIM3'
-      IF( ( nn_ice == 3 ) .AND. ( lk_cpl ) .AND. ( ( nn_limflx == -1 ) .OR. ( nn_limflx == 1 ) ) )   &
+      IF( ( nn_ice == 3 ) .AND. ( ln_cpl ) .AND. ( ( nn_limflx == -1 ) .OR. ( nn_limflx == 1 ) ) )   &
          &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
-      IF( ( nn_ice == 3 ) .AND. ( .NOT. lk_cpl ) .AND. ( nn_limflx == 2 ) )   &
+      IF( ( nn_ice == 3 ) .AND. ( .NOT. ln_cpl ) .AND. ( nn_limflx == 2 ) )   &
          &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in forced mode cannot be 2' )
 
       IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
 
-      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) )   &
+      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   &
          &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
       
-      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
-         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
-
       IF ( ln_wave ) THEN
       !Activated wave module but neither drag nor stokes drift activated
@@ -227 +236 @@
       ENDIF 
       !                          ! Choice of the Surface Boudary Condition (set nsbc)
+      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
+      !
       icpt = 0
-      IF( ln_ana          ) THEN   ;   nsbc = jp_ana    ; icpt = icpt + 1   ;   ENDIF       ! analytical      formulation
-      IF( ln_flx          ) THEN   ;   nsbc = jp_flx    ; icpt = icpt + 1   ;   ENDIF       ! flux            formulation
-      IF( ln_blk_clio     ) THEN   ;   nsbc = jp_clio   ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk       formulation
-      IF( ln_blk_core     ) THEN   ;   nsbc = jp_core   ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk       formulation
-      IF( ln_blk_mfs      ) THEN   ;   nsbc = jp_mfs    ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk       formulation
-      IF( lk_cpl          ) THEN   ;   nsbc = jp_cpl    ; icpt = icpt + 1   ;   ENDIF       ! Coupled         formulation
-      IF( cp_cfg == 'gyre') THEN   ;   nsbc = jp_gyre                       ;   ENDIF       ! GYRE analytical formulation
-      IF( lk_esopa        )            nsbc = jp_esopa                                      ! esopa test, ALL formulations
+      IF( ln_ana          ) THEN   ;   nsbc = jp_ana     ; icpt = icpt + 1   ;   ENDIF       ! analytical           formulation
+      IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
+      IF( ln_blk_clio     ) THEN   ;   nsbc = jp_clio    ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk            formulation
+      IF( ln_blk_core     ) THEN   ;   nsbc = jp_core    ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk            formulation
+      IF( ln_blk_mfs      ) THEN   ;   nsbc = jp_mfs     ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk            formulation
+      IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
+      IF( cp_cfg == 'gyre') THEN   ;   nsbc = jp_gyre                        ;   ENDIF       ! GYRE analytical      formulation
+      IF( nn_components == jp_iam_opa )   &
+         &                  THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
+      IF( lk_esopa        )            nsbc = jp_esopa                                       ! esopa test, ALL formulations
       !
       IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN
@@ -246 +259 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nsbc == jp_esopa )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
-         IF( nsbc == jp_gyre  )   WRITE(numout,*) '              GYRE analytical formulation'
-         IF( nsbc == jp_ana   )   WRITE(numout,*) '              analytical formulation'
-         IF( nsbc == jp_flx   )   WRITE(numout,*) '              flux formulation'
-         IF( nsbc == jp_clio  )   WRITE(numout,*) '              CLIO bulk formulation'
-         IF( nsbc == jp_core  )   WRITE(numout,*) '              CORE bulk formulation'
-         IF( nsbc == jp_cpl   )   WRITE(numout,*) '              coupled formulation'
-         IF( nsbc == jp_mfs   )   WRITE(numout,*) '              MFS Bulk formulation'
-      ENDIF
-      !
+         IF( nsbc == jp_esopa   )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
+         IF( nsbc == jp_gyre    )   WRITE(numout,*) '              GYRE analytical formulation'
+         IF( nsbc == jp_ana     )   WRITE(numout,*) '              analytical formulation'
+         IF( nsbc == jp_flx     )   WRITE(numout,*) '              flux formulation'
+         IF( nsbc == jp_clio    )   WRITE(numout,*) '              CLIO bulk formulation'
+         IF( nsbc == jp_core    )   WRITE(numout,*) '              CORE bulk formulation'
+         IF( nsbc == jp_purecpl )   WRITE(numout,*) '              pure coupled formulation'
+         IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
+         IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
+         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation'
+         IF( nn_components/= jp_iam_nemo )  &
+            &                       WRITE(numout,*) '              + OASIS coupled SAS'
+      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
+
+!     nn_fsbc initialization if OPA-SAS coupling via OASIS
+!     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
+      IF ( nn_components /= jp_iam_nemo ) THEN
+
+         IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
+         IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt)
+         !
+         IF(lwp)THEN
+            WRITE(numout,*)
+            WRITE(numout,*)"   OPA-SAS coupled via OASIS : nn_fsbc re-defined from OASIS namcouple ", nn_fsbc
+            WRITE(numout,*)
+         ENDIF
+      ENDIF
+
+      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
+          MOD( nstock             , nn_fsbc) /= 0 ) THEN 
+         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
+            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
+         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
+      ENDIF
+      !
+      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
+         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
+      !
+      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
+         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
+
                                CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
       !
@@ -265 +312 @@
 
       IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
-      !
-      IF( nsbc   == jp_cpl )   CALL sbc_cpl_init (nn_ice)      ! OASIS initialisation. must be done before first time step
       
    END SUBROUTINE sbc_init
@@ -310 +355 @@
                                                          ! (caution called before sbc_ssm)
       !
-      CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
-      !                                                  ! averaged over nf_sbc time-step
+      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      !                                                        ! averaged over nf_sbc time-step
 
       IF (ln_wave) CALL sbc_wave( kt )
@@ -322 +367 @@
       CASE( jp_flx   )   ;   CALL sbc_flx     ( kt )                    ! flux formulation
       CASE( jp_clio  )   ;   CALL sbc_blk_clio( kt )                    ! bulk formulation : CLIO for the ocean
-      CASE( jp_core  )   ;   CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
-      CASE( jp_cpl   )   ;   CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! coupled formulation
+      CASE( jp_core  )   
+         IF( nn_components == jp_iam_sas ) &
+            &                CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA 
+                             CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
+                                                                        ! from oce: sea surface variables (sst_m, sss_m,  ssu_m,  ssv_m)
+      CASE( jp_purecpl )  ;  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation
+                                                                        !
       CASE( jp_mfs   )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean
+      CASE( jp_none  ) 
+         IF( nn_components == jp_iam_opa ) &
+                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
       CASE( jp_esopa )                                
                              CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
@@ -334 +387 @@
       END SELECT
 
+      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
+
+
       !                                            !==  Misc. Options  ==!
       
@@ -356 +412 @@
       !                                                           ! (update freshwater fluxes)
 !RBbug do not understand why see ticket 667
-      !clem-bugsal CALL lbc_lnk( emp, 'T', 1. )
+!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
+      CALL lbc_lnk( emp, 'T', 1. )
       !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
@@ -397 +454 @@
          ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
          CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
-         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b' , sfx )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  )
       ENDIF
 
@@ -412 +469 @@
          CALL iom_put( "qns"   , qns        )                   ! solar heat flux
          CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
-         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
+         IF( nn_ice > 0 .OR. nn_components == jp_iam_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
          CALL iom_put( "taum"  , taum       )                   ! wind stress module 
          CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

@@ -32 +32 @@
 
    PUBLIC   sbc_rnf       ! routine call in sbcmod module
-   PUBLIC   sbc_rnf_div   ! routine called in sshwzv module
+   PUBLIC   sbc_rnf_div   ! routine called in divcurl module
    PUBLIC   sbc_rnf_alloc ! routine call in sbcmod module
    PUBLIC   sbc_rnf_init  ! (PUBLIC for TAM)
@@ -44 +44 @@
    LOGICAL           , PUBLIC ::   ln_rnf_tem      !: temperature river runoffs attribute specified in a file
    LOGICAL           , PUBLIC ::   ln_rnf_sal      !: salinity    river runoffs attribute specified in a file
-   LOGICAL           , PUBLIC ::   ln_rnf_emp      !: runoffs into a file to be read or already into precipitation
    TYPE(FLD_N)       , PUBLIC ::   sn_rnf          !: information about the runoff file to be read
    TYPE(FLD_N)       , PUBLIC ::   sn_cnf          !: information about the runoff mouth file to be read
@@ -54 +53 @@
    REAL(wp)          , PUBLIC ::   rn_avt_rnf      !: runoffs, value of the additional vertical mixing coef. [m2/s]
    REAL(wp)          , PUBLIC ::   rn_rfact        !: multiplicative factor for runoff
+
+   LOGICAL           , PUBLIC ::   l_rnfcpl = .false.       ! runoffs recieved from oasis
 
    INTEGER , PUBLIC  ::   nkrnf = 0         !: nb of levels over which Kz is increased at river mouths
@@ -117 +118 @@
       ENDIF
 
-      !                                                   !-------------------!
-      IF( .NOT. ln_rnf_emp ) THEN                         !   Update runoff   !
-         !                                                !-------------------!
-         !
-                             CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt
-         IF( ln_rnf_tem  )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required
-         IF( ln_rnf_sal  )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
-         !
-         ! 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 )   THEN
-            WHERE( 40._wp < gphit(:,:) .AND. gphit(:,:) < 65._wp )
-               sf_rnf(1)%fnow(:,:,1) = 0.85 * sf_rnf(1)%fnow(:,:,1)
+      !                                            !-------------------!
+      !                                            !   Update runoff   !
+      !                                            !-------------------!
+      !
+      IF( .NOT. l_rnfcpl )   CALL fld_read ( kt, nn_fsbc, sf_rnf   )    ! Read Runoffs data and provide it at kt
+      IF(   ln_rnf_tem   )   CALL fld_read ( kt, nn_fsbc, sf_t_rnf )    ! idem for runoffs temperature if required
+      IF(   ln_rnf_sal   )   CALL fld_read ( kt, nn_fsbc, sf_s_rnf )    ! idem for runoffs salinity    if required
+      !
+      ! 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
+         !
+         !                                                     ! set temperature & salinity content of runoffs
+         IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
+            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
+               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
             END WHERE
-         ENDIF
-         !
-         IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
-            !
-            rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt
-            !
-            !                                                     ! set temperature & salinity content of runoffs
-            IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data
-               rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
-               WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature
-                   rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
-               END WHERE
-               WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-                   ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study)
-                   rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp
-               END WHERE
-            ELSE                                                        ! use SST as runoffs temperature
-               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
-            ENDIF
-            !                                                           ! use runoffs salinity data
-            IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
-            !                                                           ! else use S=0 for runoffs (done one for all in the init)
-            IF ( ANY( rnf(:,:) < 0._wp ) ) z_err=1
-            IF(lk_mpp) CALL mpp_sum(z_err)
-            IF( z_err > 0 ) CALL ctl_stop( 'sbc_rnf : negative runnoff values exist' )
-            !
-            CALL iom_put( "runoffs", rnf )         ! output runoffs arrays
-         ENDIF
-         !
-      ENDIF
-      !
+            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
+               ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study)
+               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp
+            END WHERE
+         ELSE                                                        ! use SST as runoffs temperature
+            rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
+         ENDIF
+         !                                                           ! use runoffs salinity data
+         IF( ln_rnf_sal )   rnf_tsc(:,:,jp_sal) = ( sf_s_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0
+         !                                                           ! else use S=0 for runoffs (done one for all in the init)
+         CALL iom_put( "runoffs", rnf )         ! output runoffs arrays
+      ENDIF
+      !
+      !                                                ! ---------------------------------------- !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
          !                                             ! ---------------------------------------- !
@@ -172 +168 @@
          ELSE                                                   !* no restart: set from nit000 values
             IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000'
-             rnf_b    (:,:  ) = rnf    (:,:  )
-             rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+            rnf_b    (:,:  ) = rnf    (:,:  )
+            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
          ENDIF
       ENDIF
@@ -187 +183 @@
          CALL iom_rstput( kt, nitrst, numrow, 'rnf_sc_b', rnf_tsc(:,:,jp_sal) )
       ENDIF
+      !
       CALL wrk_dealloc( jpi,jpj, ztfrz)
       !
@@ -265 +262 @@
       REAL(wp), DIMENSION(:,:  ), ALLOCATABLE :: zrnf
       !
-      NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   &
+      NAMELIST/namsbc_rnf/ cn_dir            , ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   &
          &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &
          &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact,     &
@@ -290 +287 @@
          WRITE(numout,*) '~~~~~~~ '
          WRITE(numout,*) '   Namelist namsbc_rnf'
-         WRITE(numout,*) '      runoff in a file to be read                ln_rnf_emp   = ', ln_rnf_emp
          WRITE(numout,*) '      specific river mouths treatment            ln_rnf_mouth = ', ln_rnf_mouth
          WRITE(numout,*) '      river mouth additional Kz                  rn_avt_rnf   = ', rn_avt_rnf
@@ -296 +292 @@
          WRITE(numout,*) '      multiplicative factor for runoff           rn_rfact     = ', rn_rfact
       ENDIF
-      !
       !                                   ! ==================
       !                                   !   Type of runoff
@@ -303 +298 @@
       IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' )
       !
-      IF( ln_rnf_emp ) THEN                     !==  runoffs directly provided in the precipitations  ==!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          runoffs directly provided in the precipitations'
-         IF( ln_rnf_depth .OR. ln_rnf_tem .OR. ln_rnf_sal .OR. ln_rnf_depth_ini ) THEN
-           CALL ctl_warn( 'runoffs already included in precipitations, so runoff (T,S, depth) attributes will not be used' )
-           ln_rnf_depth = .FALSE.   ;   ln_rnf_tem = .FALSE.   ;   ln_rnf_sal = .FALSE.  ;   ln_rnf_depth_ini = .FALSE.
-         ENDIF
-         !
-      ELSE                                      !==  runoffs read in a file : set sf_rnf structure  ==!
-         !
+      IF( .NOT. l_rnfcpl ) THEN                    
          ALLOCATE( sf_rnf(1), STAT=ierror )         ! Create sf_rnf structure (runoff inflow)
          IF(lwp) WRITE(numout,*)
@@ -321 +307 @@
          ALLOCATE( sf_rnf(1)%fnow(jpi,jpj,1)   )
          IF( sn_rnf%ln_tint ) ALLOCATE( sf_rnf(1)%fdta(jpi,jpj,1,2) )
-         !                                          ! fill sf_rnf with the namelist (sn_rnf) and control print
          CALL fld_fill( sf_rnf, (/ sn_rnf /), cn_dir, 'sbc_rnf_init', 'read runoffs data', 'namsbc_rnf' )
-         !
-         IF( ln_rnf_tem ) THEN                      ! Create (if required) sf_t_rnf structure
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '          runoffs temperatures read in a file'
-            ALLOCATE( sf_t_rnf(1), STAT=ierror  )
-            IF( ierror > 0 ) THEN
-               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_t_rnf structure' )   ;   RETURN
-            ENDIF
-            ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj,1)   )
-            IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )
-         ENDIF
-         !
-         IF( ln_rnf_sal  ) THEN                     ! Create (if required) sf_s_rnf and sf_t_rnf structures
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '          runoffs salinities read in a file'
-            ALLOCATE( sf_s_rnf(1), STAT=ierror  )
-            IF( ierror > 0 ) THEN
-               CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_s_rnf structure' )   ;   RETURN
-            ENDIF
-            ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1)   )
-            IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) )
-            CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )
-         ENDIF
-         !
-         IF( ln_rnf_depth ) THEN                    ! depth of runoffs set from a file
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '          runoffs depth read in a file'
-            rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname )
-            IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year 
-               IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month 
-            ENDIF 
-            CALL iom_open ( rn_dep_file, inum )                           ! open file
-            CALL iom_get  ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array
-            CALL iom_close( inum )                                        ! close file
-            !
-            nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( h_rnf(ji,jj) > 0._wp ) THEN
-                     jk = 2
-                     DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 ;  END DO
-                     nk_rnf(ji,jj) = jk
-                  ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
-                  ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
-                  ELSE
-                     CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
-                     WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
-                  ENDIF
+      ENDIF
+      !
+      IF( ln_rnf_tem ) THEN                      ! Create (if required) sf_t_rnf structure
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs temperatures read in a file'
+         ALLOCATE( sf_t_rnf(1), STAT=ierror  )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_t_rnf structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_t_rnf(1)%fnow(jpi,jpj,1)   )
+         IF( sn_t_rnf%ln_tint ) ALLOCATE( sf_t_rnf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill (sf_t_rnf, (/ sn_t_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff temperature data', 'namsbc_rnf' )
+      ENDIF
+      !
+      IF( ln_rnf_sal  ) THEN                     ! Create (if required) sf_s_rnf and sf_t_rnf structures
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs salinities read in a file'
+         ALLOCATE( sf_s_rnf(1), STAT=ierror  )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'sbc_rnf_init: unable to allocate sf_s_rnf structure' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_s_rnf(1)%fnow(jpi,jpj,1)   )
+         IF( sn_s_rnf%ln_tint ) ALLOCATE( sf_s_rnf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill (sf_s_rnf, (/ sn_s_rnf /), cn_dir, 'sbc_rnf_init', 'read runoff salinity data', 'namsbc_rnf' )
+      ENDIF
+      !
+      IF( ln_rnf_depth ) THEN                    ! depth of runoffs set from a file
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          runoffs depth read in a file'
+         rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname )
+         IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year 
+            IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month 
+         ENDIF
+         CALL iom_open ( rn_dep_file, inum )                           ! open file
+         CALL iom_get  ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array
+         CALL iom_close( inum )                                        ! close file
+         !
+         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( h_rnf(ji,jj) > 0._wp ) THEN
+                  jk = 2
+                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  END DO
+                  nk_rnf(ji,jj) = jk
+               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
+               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
+               ELSE
+                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
+                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
+               ENDIF
+            END DO
+         END DO
+         DO jj = 1, jpj                                ! set the associated depth
+            DO ji = 1, jpi
+               h_rnf(ji,jj) = 0._wp
+               DO jk = 1, nk_rnf(ji,jj)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
                END DO
             END DO
-            DO jj = 1, jpj                                ! set the associated depth
-               DO ji = 1, jpi
-                  h_rnf(ji,jj) = 0._wp
-                  DO jk = 1, nk_rnf(ji,jj)
-                     h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
+         END DO
+         !
+      ELSE IF( ln_rnf_depth_ini ) THEN           ! runoffs applied at the surface
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '    depth of runoff computed once from max value of runoff'
+         IF(lwp) WRITE(numout,*) '    max value of the runoff climatologie (over global domain) rn_rnf_max = ', rn_rnf_max
+         IF(lwp) WRITE(numout,*) '    depth over which runoffs is spread                        rn_dep_max = ', rn_dep_max
+         IF(lwp) WRITE(numout,*) '     create (=1) a runoff depth file or not (=0)      nn_rnf_depth_file  = ', nn_rnf_depth_file
+
+         CALL iom_open( TRIM( sn_rnf%clname ), inum )    !  open runoff file
+         CALL iom_gettime( inum, zrec, kntime=nbrec)
+         ALLOCATE( zrnfcl(jpi,jpj,nbrec) )     ;      ALLOCATE( zrnf(jpi,jpj) )
+         DO jm = 1, nbrec
+            CALL iom_get( inum, jpdom_data, TRIM( sn_rnf%clvar ), zrnfcl(:,:,jm), jm )
+         END DO
+         CALL iom_close( inum )
+         zrnf(:,:) = MAXVAL( zrnfcl(:,:,:), DIM=3 )   !  maximum value in time
+         DEALLOCATE( zrnfcl )
+         !
+         h_rnf(:,:) = 1.
+         !
+         zacoef = rn_dep_max / rn_rnf_max            ! coef of linear relation between runoff and its depth (150m for max of runoff)
+         !
+         WHERE( zrnf(:,:) > 0._wp )  h_rnf(:,:) = zacoef * zrnf(:,:)   ! compute depth for all runoffs
+         !
+         DO jj = 1, jpj                     ! take in account min depth of ocean rn_hmin
+            DO ji = 1, jpi
+               IF( zrnf(ji,jj) > 0._wp ) THEN
+                  jk = mbkt(ji,jj)
+                  h_rnf(ji,jj) = MIN( h_rnf(ji,jj), gdept_0(ji,jj,jk ) )
+               ENDIF
+            END DO
+         END DO
+         !
+         nk_rnf(:,:) = 0                       ! number of levels on which runoffs are distributed
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zrnf(ji,jj) > 0._wp ) THEN
+                  jk = 2
+                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
                   END DO
+                  nk_rnf(ji,jj) = jk
+               ELSE
+                  nk_rnf(ji,jj) = 1
+               ENDIF
+            END DO
+         END DO
+         !
+         DEALLOCATE( zrnf )
+         !
+         DO jj = 1, jpj                                ! set the associated depth
+            DO ji = 1, jpi
+               h_rnf(ji,jj) = 0._wp
+               DO jk = 1, nk_rnf(ji,jj)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
                END DO
             END DO
-            !
-         ELSE IF( ln_rnf_depth_ini ) THEN           ! runoffs applied at the surface
-            !
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '    depth of runoff computed once from max value of runoff'
-            IF(lwp) WRITE(numout,*) '    max value of the runoff climatologie (over global domain) rn_rnf_max = ', rn_rnf_max
-            IF(lwp) WRITE(numout,*) '    depth over which runoffs is spread                        rn_dep_max = ', rn_dep_max
-            IF(lwp) WRITE(numout,*) '     create (=1) a runoff depth file or not (=0)      nn_rnf_depth_file  = ', nn_rnf_depth_file
-
-            CALL iom_open( TRIM( sn_rnf%clname ), inum )    !  open runoff file
-            CALL iom_gettime( inum, zrec, kntime=nbrec)
-            ALLOCATE( zrnfcl(jpi,jpj,nbrec) )     ;      ALLOCATE( zrnf(jpi,jpj) )
-            DO jm = 1, nbrec
-               CALL iom_get( inum, jpdom_data, TRIM( sn_rnf%clvar ), zrnfcl(:,:,jm), jm )
-            END DO
-            CALL iom_close( inum )
-            zrnf(:,:) = MAXVAL( zrnfcl(:,:,:), DIM=3 )   !  maximum value in time
-            DEALLOCATE( zrnfcl )
-            !
-            h_rnf(:,:) = 1.
-            !
-            zacoef = rn_dep_max / rn_rnf_max            ! coef of linear relation between runoff and its depth (150m for max of runoff)
-            !
-            WHERE( zrnf(:,:) > 0._wp )  h_rnf(:,:) = zacoef * zrnf(:,:)   ! compute depth for all runoffs
-            !
-            DO jj = 1, jpj                     ! take in account min depth of ocean rn_hmin
-               DO ji = 1, jpi
-                  IF( zrnf(ji,jj) > 0._wp ) THEN
-                     jk = mbkt(ji,jj)
-                     h_rnf(ji,jj) = MIN( h_rnf(ji,jj), gdept_0(ji,jj,jk ) )
-                  ENDIF
-               END DO
-            END DO
-            !
-            nk_rnf(:,:) = 0                       ! number of levels on which runoffs are distributed
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                   IF( zrnf(ji,jj) > 0._wp ) THEN
-                     jk = 2
-                     DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1 ;  END DO
-                     nk_rnf(ji,jj) = jk
-                   ELSE
-                     nk_rnf(ji,jj) = 1
-                   ENDIF
-                END DO
-            END DO
-            !
-            DEALLOCATE( zrnf )
-            !
-            DO jj = 1, jpj                                ! set the associated depth
-               DO ji = 1, jpi
-                  h_rnf(ji,jj) = 0._wp
-                  DO jk = 1, nk_rnf(ji,jj)
-                     h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-            !
-            IF( nn_rnf_depth_file == 1 ) THEN      !  save  output nb levels for runoff
-               IF(lwp) WRITE(numout,*) '              create runoff depht file'
-               CALL iom_open  ( TRIM( sn_dep_rnf%clname ), inum, ldwrt = .TRUE., kiolib = jprstlib )
-               CALL iom_rstput( 0, 0, inum, 'rodepth', h_rnf )
-               CALL iom_close ( inum )
-            ENDIF
-         ELSE                                       ! runoffs applied at the surface
-            nk_rnf(:,:) = 1
-            h_rnf (:,:) = fse3t(:,:,1)
-         ENDIF
-         !
+         END DO
+         !
+         IF( nn_rnf_depth_file == 1 ) THEN      !  save  output nb levels for runoff
+            IF(lwp) WRITE(numout,*) '              create runoff depht file'
+            CALL iom_open  ( TRIM( sn_dep_rnf%clname ), inum, ldwrt = .TRUE., kiolib = jprstlib )
+            CALL iom_rstput( 0, 0, inum, 'rodepth', h_rnf )
+            CALL iom_close ( inum )
+         ENDIF
+      ELSE                                       ! runoffs applied at the surface
+         nk_rnf(:,:) = 1
+         h_rnf (:,:) = fse3t(:,:,1)
       ENDIF
       !
@@ -471 +457 @@
          IF( rn_hrnf > 0._wp ) THEN
             nkrnf = 2
-            DO WHILE( nkrnf /= jpkm1 .AND. gdepw_1d(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1   ;   END DO
+            DO WHILE( nkrnf /= jpkm1 .AND. gdepw_1d(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1
+            END DO
             IF( ln_sco )   CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' )
          ENDIF

trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

@@ -58 +58 @@
       REAL(wp) ::   zcoef, zf_sbc       ! local scalar
       REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts
-      REAL(wp), DIMENSION(jpi,jpj)      :: zub, zvb,zdep
       !!---------------------------------------------------------------------
-      
+
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
       DO jj = 1, jpj
@@ -68 +67 @@
          END DO
       END DO
-      zub(:,:)        = ub (:,:,1       )
-      zvb(:,:)        = vb (:,:,1       )
-      !
-      IF( lk_vvl ) THEN
-         zdep(:,:) = fse3t_n(:,:,1)
-      ENDIF
-      !                                                   ! ---------------------------------------- !
+      !
       IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        !
          !                                                ! ---------------------------------------- !
-         ssu_m(:,:) = zub(:,:)
-         ssv_m(:,:) = zvb(:,:)
+         ssu_m(:,:) = ub(:,:,1)
+         ssv_m(:,:) = vb(:,:,1)
          IF( ln_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
          ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem)
@@ -88 +81 @@
          ENDIF
          !
-         IF( lk_vvl )   fse3t_m(:,:) = zdep(:,:)
+         IF( lk_vvl )   e3t_m(:,:) = fse3t_n(:,:,1)
+         !
+         frq_m(:,:) = fraqsr_1lev(:,:)
          !
       ELSE
@@ -97 +92 @@
             IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields initialised to instantaneous values'
             zcoef = REAL( nn_fsbc - 1, wp )
-            ssu_m(:,:) = zcoef * zub(:,:)
-            ssv_m(:,:) = zcoef * zvb(:,:)
+            ssu_m(:,:) = zcoef * ub(:,:,1)
+            ssv_m(:,:) = zcoef * vb(:,:,1)
             IF( ln_useCT )  THEN    ;   sst_m(:,:) = zcoef * eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
             ELSE                    ;   sst_m(:,:) = zcoef * zts(:,:,jp_tem)
@@ -108 +103 @@
             ENDIF
             !
-            IF( lk_vvl )   fse3t_m(:,:) = zcoef * zdep(:,:)
+            IF( lk_vvl )   e3t_m(:,:) = zcoef * fse3t_n(:,:,1)
+            !
+            frq_m(:,:) = zcoef * fraqsr_1lev(:,:)
             !                                             ! ---------------------------------------- !
          ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   !
@@ -117 +114 @@
             sss_m(:,:) = 0.e0
             ssh_m(:,:) = 0.e0
-            IF( lk_vvl )   fse3t_m(:,:) = 0.e0
+            IF( lk_vvl )   e3t_m(:,:) = 0.e0
+            frq_m(:,:) = 0.e0
          ENDIF
          !                                                ! ---------------------------------------- !
          !                                                !        Cumulate at each time step        !
          !                                                ! ---------------------------------------- !
-         ssu_m(:,:) = ssu_m(:,:) + zub(:,:)
-         ssv_m(:,:) = ssv_m(:,:) + zvb(:,:)
+         ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1)
+         ssv_m(:,:) = ssv_m(:,:) + vb(:,:,1)
          IF( ln_useCT )  THEN    ;   sst_m(:,:) = sst_m(:,:) + eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
          ELSE                    ;   sst_m(:,:) = sst_m(:,:) + zts(:,:,jp_tem)
@@ -133 +131 @@
          ENDIF
          !
-         IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) + zdep(:,:)
+         IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1)
+         !
+         frq_m(:,:) =   frq_m(:,:) + fraqsr_1lev(:,:)
 
          !                                                ! ---------------------------------------- !
@@ -144 +144 @@
             ssv_m(:,:) = ssv_m(:,:) * zcoef           !
             ssh_m(:,:) = ssh_m(:,:) * zcoef           ! mean SSH             [m]
-            IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m]
+            IF( lk_vvl )   e3t_m(:,:) = fse3t_m(:,:) * zcoef   ! mean vertical scale factor [m]
+            frq_m(:,:) = frq_m(:,:) * zcoef   ! mean fraction of solar net radiation absorbed in the 1st T level [-]
             !
          ENDIF
@@ -161 +162 @@
             CALL iom_rstput( kt, nitrst, numrow, 'sss_m'  , sss_m  )
             CALL iom_rstput( kt, nitrst, numrow, 'ssh_m'  , ssh_m  )
-            IF( lk_vvl ) THEN
-               CALL iom_rstput( kt, nitrst, numrow, 'fse3t_m'  , fse3t_m(:,:)  )
-            END IF
-            !
-         ENDIF
-         !
+            IF( lk_vvl )   CALL iom_rstput( kt, nitrst, numrow, 'e3t_m'  , e3t_m  )
+            CALL iom_rstput( kt, nitrst, numrow, 'frq_m'  , frq_m  )
+            !
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( MOD( kt - 1 , nn_fsbc ) == 0 ) THEN          !   Mean value at each nn_fsbc time-step   !
+         CALL iom_put( 'ssu_m', ssu_m )
+         CALL iom_put( 'ssv_m', ssv_m )
+         CALL iom_put( 'sst_m', sst_m )
+         CALL iom_put( 'sss_m', sss_m )
+         CALL iom_put( 'ssh_m', ssh_m )
+         IF( lk_vvl )   CALL iom_put( 'e3t_m', e3t_m )
+         CALL iom_put( 'frq_m', frq_m )
       ENDIF
       !
@@ -202 +212 @@
             CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m  )   !   "         "    salinity    (T-point)
             CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m  )   !   "         "    height      (T-point)
-            IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'fse3t_m', fse3t_m(:,:) )
+            IF( lk_vvl ) CALL iom_get( numror, jpdom_autoglo, 'e3t_m', e3t_m )
+            ! fraction of solar net radiation absorbed in 1st T level
+            IF( iom_varid( numror, 'frq_m', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'frq_m'  , frq_m  )
+            ELSE
+               frq_m(:,:) = 1._wp   ! default definition
+            ENDIF
             !
             IF( zf_sbc /= REAL( nn_fsbc, wp ) ) THEN      ! nn_fsbc has changed between 2 runs
@@ -213 +229 @@
                sss_m(:,:) = zcoef * sss_m(:,:)
                ssh_m(:,:) = zcoef * ssh_m(:,:)
-               IF( lk_vvl ) fse3t_m(:,:) = zcoef * fse3t_m(:,:)
+               IF( lk_vvl )   e3t_m(:,:) = zcoef * fse3t_m(:,:)
+               frq_m(:,:) = zcoef * frq_m(:,:)
             ELSE
                IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields read in the ocean restart file'
@@ -220 +237 @@
       ENDIF
       !
+      IF( .NOT. l_ssm_mean ) THEN   ! default initialisation. needed by lim_istate
+         !
+         IF(lwp) WRITE(numout,*) '          default initialisation of ss?_m arrays'
+         ssu_m(:,:) = ub(:,:,1)
+         ssv_m(:,:) = vb(:,:,1)
+         IF( ln_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
+         ELSE                    ;   sst_m(:,:) = tsn(:,:,1,jp_tem)
+         ENDIF
+         sss_m(:,:) = tsn(:,:,1,jp_sal)
+         ssh_m(:,:) = sshn(:,:)
+         IF( lk_vvl )   e3t_m(:,:) = fse3t_n(:,:,1)
+         frq_m(:,:) = 1._wp
+         !
+      ENDIF
+      !
    END SUBROUTINE sbc_ssm_init
 

trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

@@ -74 +74 @@
    PUBLIC   eos_init       ! called by istate module
 
-   !                                          !!* Namelist (nameos) *
-   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
-   LOGICAL , PUBLIC ::   ln_useCT  = .FALSE.  ! determine if eos_pt_from_ct is used to compute sst_m
+   !                                !!* Namelist (nameos) *
+   INTEGER , PUBLIC ::   nn_eos     ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
+   LOGICAL , PUBLIC ::   ln_useCT   ! determine if eos_pt_from_ct is used to compute sst_m
 
    !                                   !!!  simplified eos coefficients
@@ -1252 +1252 @@
             WRITE(numout,*) '             model uses Conservative Temperature'
             WRITE(numout,*) '             Important: model must be initialized with CT and SA fields'
+         ELSE
+            WRITE(numout,*) '             model does not use Conservative Temperature'
          ENDIF
       ENDIF

trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -32 +32 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-   USE sbc_ice, ONLY : lk_lim3
 
    IMPLICIT NONE
@@ -38 +37 @@
 
    PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T)
-   PUBLIC   tra_qsr_init  ! routine called by opa.F90
+   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90
 
    !                                 !!* Namelist namtra_qsr: penetrative solar radiation
@@ -50 +49 @@
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
-   
+ 
    ! Module variables
    REAL(wp) ::   xsi0r                           !: inverse of rn_si0
@@ -165 +164 @@
          CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution
          ! clem: store attenuation coefficient of the first ocean level
-         IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
+         IF ( ln_qsr_ice ) THEN
             DO jj = 1, jpj
                DO ji = 1, jpi
                   IF ( qsr(ji,jj) /= 0._wp ) THEN
                      fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) )
+                  ELSE
+                     fraqsr_1lev(ji,jj) = 1.
                   ENDIF
                END DO
@@ -233 +234 @@
                END DO
                ! clem: store attenuation coefficient of the first ocean level
-               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
+               IF ( ln_qsr_ice ) THEN
                   DO jj = 1, jpj
                      DO ji = 1, jpi
@@ -256 +257 @@
                END DO
                ! clem: store attenuation coefficient of the first ocean level
-               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
+               IF ( ln_qsr_ice ) THEN
                   fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
                ENDIF
@@ -279 +280 @@
                END DO
                ! clem: store attenuation coefficient of the first ocean level
-               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
+               IF ( ln_qsr_ice ) THEN
                   DO jj = 1, jpj
                      DO ji = 1, jpi
@@ -298 +299 @@
                END DO
                ! clem: store attenuation coefficient of the first ocean level
-               IF ( lk_lim3 .AND. ln_qsr_ice ) THEN
+               IF ( ln_qsr_ice ) THEN
                   fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp
                ENDIF
@@ -324 +325 @@
             &                    'at it= ', kt,' date= ', ndastp
          IF(lwp) WRITE(numout,*) '~~~~'
-         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc )
+         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      )
+         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm 
          !
       ENDIF
@@ -379 +381 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init')
-      !
-      ! Default value for fraqsr_1lev
-      IF( .NOT. ln_rstart ) THEN
-         fraqsr_1lev(:,:) = 1._wp
-      ENDIF
       !
       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
@@ -412 +409 @@
          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0
          WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1
-         WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice    
       ENDIF
 
@@ -564 +560 @@
       ENDIF
       !
+      ! initialisation of fraqsr_1lev used in sbcssm
+      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN
+         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  )
+      ELSE
+         fraqsr_1lev(:,:) = 1._wp   ! default definition
+      ENDIF
+      !
       CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        ) 
       CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea ) 

trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -761 +761 @@
       IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
       IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
-      IF( nn_etau == 3 .AND. .NOT. lk_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
+      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
 
       IF( ln_mxl0 ) THEN

trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -82 +82 @@
    USE crsini          ! initialise grid coarsening utility
    USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 
+   USE sbc_oce, ONLY: lk_oasis
    USE stopar
    USE stopts
@@ -197 +198 @@
 #if defined key_iomput
       CALL xios_finalize                ! end mpp communications with xios
-      IF( lk_cpl ) CALL cpl_finalize    ! end coupling and mpp communications with OASIS
+      IF( lk_oasis ) CALL cpl_finalize    ! end coupling and mpp communications with OASIS
 #else
-      IF( lk_cpl ) THEN 
+      IF( lk_oasis ) THEN 
          CALL cpl_finalize              ! end coupling and mpp communications with OASIS
       ELSE
@@ -228 +229 @@
       !
       cltxt = ''
+      cxios_context = 'nemo'
       !
       !                             ! Open reference namelist and configuration namelist files
@@ -274 +276 @@
 #if defined key_iomput
       IF( Agrif_Root() ) THEN
-         IF( lk_cpl ) THEN
-            CALL cpl_init( ilocal_comm )                               ! nemo local communicator given by oasis
-            CALL xios_initialize( "oceanx",local_comm=ilocal_comm )    ! send nemo communicator to xios
+         IF( lk_oasis ) THEN
+            CALL cpl_init( "oceanx", ilocal_comm )                     ! nemo local communicator given by oasis
+            CALL xios_initialize( "not used",local_comm=ilocal_comm )    ! send nemo communicator to xios
          ELSE
-            CALL  xios_initialize( "nemo",return_comm=ilocal_comm )    ! nemo local communicator given by xios
+            CALL  xios_initialize( "for_xios_mpi_id",return_comm=ilocal_comm )    ! nemo local communicator given by xios
          ENDIF
       ENDIF
-      narea = mynode( cltxt, numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )   ! Nodes selection
+      ! Nodes selection (control print return in cltxt)
+      narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
 #else
-      IF( lk_cpl ) THEN
+      IF( lk_oasis ) THEN
          IF( Agrif_Root() ) THEN
-            CALL cpl_init( ilocal_comm )                       ! nemo local communicator given by oasis
+            CALL cpl_init( "oceanx", ilocal_comm )                      ! nemo local communicator given by oasis
          ENDIF
-         narea = mynode( cltxt, numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )   ! Nodes selection (control print return in cltxt)
+         ! Nodes selection (control print return in cltxt)
+         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop, ilocal_comm )
       ELSE
          ilocal_comm = 0
-         narea = mynode( cltxt, numnam_ref, numnam_cfg, numond , nstop )                ! Nodes selection (control print return in cltxt)
+         ! Nodes selection (control print return in cltxt)
+         narea = mynode( cltxt, 'output.namelist.dyn', numnam_ref, numnam_cfg, numond , nstop )
       ENDIF
 #endif

trunk/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -83 +83 @@
       IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
 # if defined key_iomput
-      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( "nemo" )
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
 # endif
 #endif
                              indic = 0           ! reset to no error condition
       IF( kstp == nit000 ) THEN
-                      CALL iom_init( "nemo" )      ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
-         IF( ln_crs ) CALL iom_init( "nemo_crs" )  ! initialize context for coarse grid
+         ! must be done after nemo_init for AGRIF+XIOS+OASIS
+                      CALL iom_init(      cxios_context          )  ! iom_put initialization
+         IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" )  ! initialize context for coarse grid
       ENDIF
 
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
-                             CALL iom_setkt( kstp - nit000 + 1, "nemo"     )   ! say to iom that we are at time step kstp
-      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, "nemo_crs" )   ! say to iom that we are at time step kstp
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell iom we are at time step kstp
+      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell iom we are at time step kstp
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -168 +169 @@
       IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
 #endif
-#if defined key_traldf_c3d && key_traldf_smag
+#if defined key_traldf_c3d && defined key_traldf_smag
                           CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
 #  endif
-#if defined key_dynldf_c3d && key_dynldf_smag
+#if defined key_dynldf_c3d && defined key_dynldf_smag
                           CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
 #  endif
@@ -225 +226 @@
       IF( lk_floats  )      CALL flo_stp( kstp )         ! drifting Floats
       IF( lk_diahth  )      CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
-      IF( .NOT. lk_cpl )    CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
+      IF( .NOT. ln_cpl )    CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
       IF( lk_diadct  )      CALL dia_dct( kstp )         ! Transports
       IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag
@@ -355 +356 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+      IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
       IF( kstp == nitend .OR. indic < 0 ) THEN 
-                      CALL iom_context_finalize( "nemo"     ) ! needed for XIOS+AGRIF
-         IF( ln_crs ) CALL iom_context_finalize( "nemo_crs" ) ! 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
       ENDIF
 #endif