New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in NEMO/branches/UKMO/NEMO_4.0.4_MEDUSA_externals_GC5/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_MEDUSA_externals_GC5/src/OCE/SBC/sbccpl.F90 @ 15806

Last change on this file since 15806 was 15806, checked in by jpalmier, 2 years ago

fix -1

File size: 177.7 KB
Line 
1MODULE sbccpl
2   !!======================================================================
3   !!                       ***  MODULE  sbccpl  ***
4   !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
5   !!======================================================================
6   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
7   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
8   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
9   !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   namsbc_cpl      : coupled formulation namlist
14   !!   sbc_cpl_init    : initialisation of the coupled exchanges
15   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
16   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
17   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
18   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
19   !!   sbc_cpl_snd     : send     fields to the atmosphere
20   !!----------------------------------------------------------------------
21   USE dom_oce         ! ocean space and time domain
22   USE sbc_oce         ! Surface boundary condition: ocean fields
23   USE trc_oce         ! share SMS/Ocean variables
24   USE sbc_ice         ! Surface boundary condition: ice fields
25   USE sbcapr          ! Stochastic param. : ???
26   USE sbcdcy          ! surface boundary condition: diurnal cycle
27   USE sbcwave         ! surface boundary condition: waves
28   USE phycst          ! physical constants
29#if defined key_si3
30   USE ice            ! ice variables
31#endif
32   USE cpl_oasis3     ! OASIS3 coupling
33   USE geo2ocean      !
34   USE oce     , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
35   USE ocealb         !
36   USE eosbn2         !
37   USE sbcrnf  , ONLY : l_rnfcpl
38   USE cpl_rnf_1d, ONLY: nn_cpl_river, cpl_rnf_1d_init, cpl_rnf_1d_to_2d   ! Variables used in 1D river outflow
39   USE sbcisf  , ONLY : l_isfcpl
40   USE par_trc , ONLY : ln_medusa
41#if defined key_cice
42   USE ice_domain_size, only: ncat
43#endif
44#if defined key_si3
45   USE icevar         ! for CALL ice_var_snwblow
46#endif
47   !
48   USE in_out_manager ! I/O manager
49   USE iom            ! NetCDF library
50   USE lib_mpp        ! distribued memory computing library
51   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
52   USE timing         ! timing
53   USE sms_medusa, ONLY: CO2Flux_out_cpl, DMS_out_cpl, chloro_out_cpl, PCO2a_in_cpl, Dust_in_cpl
54
55#if defined key_oasis3 
56   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut 
57#endif
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   sbc_cpl_init      ! routine called by sbcmod.F90
63   PUBLIC   sbc_cpl_rcv       ! routine called by icestp.F90
64   PUBLIC   sbc_cpl_snd       ! routine called by step.F90
65   PUBLIC   sbc_cpl_ice_tau   ! routine called by icestp.F90
66   PUBLIC   sbc_cpl_ice_flx   ! routine called by icestp.F90
67   PUBLIC   sbc_cpl_alloc     ! routine called in sbcice_cice.F90
68
69   INTEGER, PARAMETER ::   jpr_otx1   =  1   ! 3 atmosphere-ocean stress components on grid 1
70   INTEGER, PARAMETER ::   jpr_oty1   =  2   !
71   INTEGER, PARAMETER ::   jpr_otz1   =  3   !
72   INTEGER, PARAMETER ::   jpr_otx2   =  4   ! 3 atmosphere-ocean stress components on grid 2
73   INTEGER, PARAMETER ::   jpr_oty2   =  5   !
74   INTEGER, PARAMETER ::   jpr_otz2   =  6   !
75   INTEGER, PARAMETER ::   jpr_itx1   =  7   ! 3 atmosphere-ice   stress components on grid 1
76   INTEGER, PARAMETER ::   jpr_ity1   =  8   !
77   INTEGER, PARAMETER ::   jpr_itz1   =  9   !
78   INTEGER, PARAMETER ::   jpr_itx2   = 10   ! 3 atmosphere-ice   stress components on grid 2
79   INTEGER, PARAMETER ::   jpr_ity2   = 11   !
80   INTEGER, PARAMETER ::   jpr_itz2   = 12   !
81   INTEGER, PARAMETER ::   jpr_qsroce = 13   ! Qsr above the ocean
82   INTEGER, PARAMETER ::   jpr_qsrice = 14   ! Qsr above the ice
83   INTEGER, PARAMETER ::   jpr_qsrmix = 15 
84   INTEGER, PARAMETER ::   jpr_qnsoce = 16   ! Qns above the ocean
85   INTEGER, PARAMETER ::   jpr_qnsice = 17   ! Qns above the ice
86   INTEGER, PARAMETER ::   jpr_qnsmix = 18
87   INTEGER, PARAMETER ::   jpr_rain   = 19   ! total liquid precipitation (rain)
88   INTEGER, PARAMETER ::   jpr_snow   = 20   ! solid precipitation over the ocean (snow)
89   INTEGER, PARAMETER ::   jpr_tevp   = 21   ! total evaporation
90   INTEGER, PARAMETER ::   jpr_ievp   = 22   ! solid evaporation (sublimation)
91   INTEGER, PARAMETER ::   jpr_sbpr   = 23   ! sublimation - liquid precipitation - solid precipitation
92   INTEGER, PARAMETER ::   jpr_semp   = 24   ! solid freshwater budget (sublimation - snow)
93   INTEGER, PARAMETER ::   jpr_oemp   = 25   ! ocean freshwater budget (evap - precip)
94   INTEGER, PARAMETER ::   jpr_w10m   = 26   ! 10m wind
95   INTEGER, PARAMETER ::   jpr_dqnsdt = 27   ! d(Q non solar)/d(temperature)
96   INTEGER, PARAMETER ::   jpr_rnf    = 28   ! runoffs
97   INTEGER, PARAMETER ::   jpr_cal    = 29   ! calving
98   INTEGER, PARAMETER ::   jpr_taum   = 30   ! wind stress module
99   INTEGER, PARAMETER ::   jpr_co2    = 31
100   INTEGER, PARAMETER ::   jpr_topm   = 32   ! topmeltn
101   INTEGER, PARAMETER ::   jpr_botm   = 33   ! botmeltn
102   INTEGER, PARAMETER ::   jpr_sflx   = 34   ! salt flux
103   INTEGER, PARAMETER ::   jpr_toce   = 35   ! ocean temperature
104   INTEGER, PARAMETER ::   jpr_soce   = 36   ! ocean salinity
105   INTEGER, PARAMETER ::   jpr_ocx1   = 37   ! ocean current on grid 1
106   INTEGER, PARAMETER ::   jpr_ocy1   = 38   !
107   INTEGER, PARAMETER ::   jpr_ssh    = 39   ! sea surface height
108   INTEGER, PARAMETER ::   jpr_fice   = 40   ! ice fraction         
109   INTEGER, PARAMETER ::   jpr_e3t1st = 41   ! first T level thickness
110   INTEGER, PARAMETER ::   jpr_fraqsr = 42   ! fraction of solar net radiation absorbed in the first ocean level
111   INTEGER, PARAMETER ::   jpr_mslp   = 43   ! mean sea level pressure
112   INTEGER, PARAMETER ::   jpr_hsig   = 44   ! Hsig
113   INTEGER, PARAMETER ::   jpr_phioc  = 45   ! Wave=>ocean energy flux
114   INTEGER, PARAMETER ::   jpr_sdrftx = 46   ! Stokes drift on grid 1
115   INTEGER, PARAMETER ::   jpr_sdrfty = 47   ! Stokes drift on grid 2
116   INTEGER, PARAMETER ::   jpr_wper   = 48   ! Mean wave period
117   INTEGER, PARAMETER ::   jpr_wnum   = 49   ! Mean wavenumber
118   INTEGER, PARAMETER ::   jpr_tauwoc = 50   ! Stress fraction adsorbed by waves
119   INTEGER, PARAMETER ::   jpr_wdrag  = 51   ! Neutral surface drag coefficient
120   INTEGER, PARAMETER ::   jpr_isf    = 52
121   INTEGER, PARAMETER ::   jpr_icb    = 53
122   INTEGER, PARAMETER ::   jpr_wfreq  = 54   ! Wave peak frequency
123   INTEGER, PARAMETER ::   jpr_tauwx  = 55   ! x component of the ocean stress from waves
124   INTEGER, PARAMETER ::   jpr_tauwy  = 56   ! y component of the ocean stress from waves
125   INTEGER, PARAMETER ::   jpr_ts_ice = 57   ! Sea ice surface temp
126   INTEGER, PARAMETER ::   jpr_grnm   = 58   ! Greenland ice mass
127   INTEGER, PARAMETER ::   jpr_antm   = 59   ! Antarctic ice mass
128   INTEGER, PARAMETER ::   jpr_rnf_1d = 60   ! 1D river runoff
129   INTEGER, PARAMETER ::   jpr_qtr    = 61   ! Transmitted solar
130   INTEGER, PARAMETER ::   jpr_atm_pco2 = 62 ! Incoming atm pCO2 flux
131   INTEGER, PARAMETER ::   jpr_atm_dust = 63 ! Incoming atm aggregate dust
132
133   INTEGER, PARAMETER ::   jprcv      = 63   ! total number of fields received
134
135   INTEGER, PARAMETER ::   jps_fice   =  1   ! ice fraction sent to the atmosphere
136   INTEGER, PARAMETER ::   jps_toce   =  2   ! ocean temperature
137   INTEGER, PARAMETER ::   jps_tice   =  3   ! ice   temperature
138   INTEGER, PARAMETER ::   jps_tmix   =  4   ! mixed temperature (ocean+ice)
139   INTEGER, PARAMETER ::   jps_albice =  5   ! ice   albedo
140   INTEGER, PARAMETER ::   jps_albmix =  6   ! mixed albedo
141   INTEGER, PARAMETER ::   jps_hice   =  7   ! ice  thickness
142   INTEGER, PARAMETER ::   jps_hsnw   =  8   ! snow thickness
143   INTEGER, PARAMETER ::   jps_ocx1   =  9   ! ocean current on grid 1
144   INTEGER, PARAMETER ::   jps_ocy1   = 10   !
145   INTEGER, PARAMETER ::   jps_ocz1   = 11   !
146   INTEGER, PARAMETER ::   jps_ivx1   = 12   ! ice   current on grid 1
147   INTEGER, PARAMETER ::   jps_ivy1   = 13   !
148   INTEGER, PARAMETER ::   jps_ivz1   = 14   !
149   INTEGER, PARAMETER ::   jps_co2    = 15
150   INTEGER, PARAMETER ::   jps_soce   = 16   ! ocean salinity
151   INTEGER, PARAMETER ::   jps_ssh    = 17   ! sea surface height
152   INTEGER, PARAMETER ::   jps_qsroce = 18   ! Qsr above the ocean
153   INTEGER, PARAMETER ::   jps_qnsoce = 19   ! Qns above the ocean
154   INTEGER, PARAMETER ::   jps_oemp   = 20   ! ocean freshwater budget (evap - precip)
155   INTEGER, PARAMETER ::   jps_sflx   = 21   ! salt flux
156   INTEGER, PARAMETER ::   jps_otx1   = 22   ! 2 atmosphere-ocean stress components on grid 1
157   INTEGER, PARAMETER ::   jps_oty1   = 23   !
158   INTEGER, PARAMETER ::   jps_rnf    = 24   ! runoffs
159   INTEGER, PARAMETER ::   jps_taum   = 25   ! wind stress module
160   INTEGER, PARAMETER ::   jps_fice2  = 26   ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
161   INTEGER, PARAMETER ::   jps_e3t1st = 27   ! first level depth (vvl)
162   INTEGER, PARAMETER ::   jps_fraqsr = 28   ! fraction of solar net radiation absorbed in the first ocean level
163   INTEGER, PARAMETER ::   jps_ficet  = 29   ! total ice fraction 
164   INTEGER, PARAMETER ::   jps_ocxw   = 30   ! currents on grid 1 
165   INTEGER, PARAMETER ::   jps_ocyw   = 31   ! currents on grid 2
166   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level
167   INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes)
168   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area fraction
169   INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness
170   INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity
171   INTEGER, PARAMETER ::   jps_sstfrz = 37   ! sea surface freezing temperature
172   INTEGER, PARAMETER ::   jps_ttilyr = 38   ! sea ice top layer temp
173   INTEGER, PARAMETER ::   jps_bio_co2 = 39    ! MEDUSA air-sea CO2 flux
174   INTEGER, PARAMETER ::   jps_bio_dms = 40    ! MEDUSA DMS surface concentration
175   INTEGER, PARAMETER ::   jps_bio_chloro = 41 ! MEDUSA chlorophyll surface concentration
176
177   INTEGER, PARAMETER ::   jpsnd      = 41   ! total number of fields sent
178
179#if ! defined key_oasis3
180   ! Dummy variables to enable compilation when oasis3 is not being used
181   INTEGER                    ::   OASIS_Sent        = -1 
182   INTEGER                    ::   OASIS_SentOut     = -1 
183   INTEGER                    ::   OASIS_ToRest      = -1 
184   INTEGER                    ::   OASIS_ToRestOut   = -1 
185#endif 
186
187   !                                  !!** namelist namsbc_cpl **
188   TYPE ::   FLD_C                     !   
189      CHARACTER(len = 32) ::   cldes      ! desciption of the coupling strategy
190      CHARACTER(len = 32) ::   clcat      ! multiple ice categories strategy
191      CHARACTER(len = 32) ::   clvref     ! reference of vector ('spherical' or 'cartesian')
192      CHARACTER(len = 32) ::   clvor      ! orientation of vector fields ('eastward-northward' or 'local grid')
193      CHARACTER(len = 32) ::   clvgrd     ! grids on which is located the vector fields
194   END TYPE FLD_C
195   !                                   ! Send to the atmosphere 
196   TYPE(FLD_C) ::   sn_snd_temp  , sn_snd_alb , sn_snd_thick, sn_snd_crt   , sn_snd_co2,  &
197      &             sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr
198   !                                   ! Received from the atmosphere
199   TYPE(FLD_C) ::   sn_snd_bio_co2, sn_snd_bio_dms, sn_snd_bio_chloro
200   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_tauw, sn_rcv_dqnsdt, sn_rcv_qsr,  &
201      &             sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf, sn_rcv_ts_ice
202   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf,      &
203                    sn_rcv_grnm, sn_rcv_antm
204   TYPE(FLD_C) ::   sn_rcv_atm_pco2, sn_rcv_atm_dust
205   ! Send to waves
206   TYPE(FLD_C) ::   sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 
207   ! Received from waves
208   TYPE(FLD_C) ::   sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, sn_rcv_tauwoc, &
209                    sn_rcv_wdrag, sn_rcv_wfreq
210   ! Transmitted solar
211   TYPE(FLD_C) ::   sn_rcv_qtr
212   !                                   ! Other namelist parameters
213   INTEGER     ::   nn_cplmodel           ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
214   LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models
215                                          !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
216   LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
217
218   TYPE ::   DYNARR     
219      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3   
220   END TYPE DYNARR
221
222   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                ! all fields recieved from the atmosphere
223
224   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky)
225
226   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]
227   REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
228
229   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:) ::   nrcvinfo           ! OASIS info argument
230
231   !! Substitution
232#  include "vectopt_loop_substitute.h90"
233   !!----------------------------------------------------------------------
234   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
235   !! $Id$
236   !! Software governed by the CeCILL license (see ./LICENSE)
237   !!----------------------------------------------------------------------
238CONTAINS
239 
240   INTEGER FUNCTION sbc_cpl_alloc()
241      !!----------------------------------------------------------------------
242      !!             ***  FUNCTION sbc_cpl_alloc  ***
243      !!----------------------------------------------------------------------
244      INTEGER :: ierr(4)
245      !!----------------------------------------------------------------------
246      ierr(:) = 0
247      !
248      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
249     
250#if ! defined key_si3 && ! defined key_cice
251      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
252#endif
253      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
254      !
255      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) ) 
256
257      sbc_cpl_alloc = MAXVAL( ierr )
258      CALL mpp_sum ( 'sbccpl', sbc_cpl_alloc )
259      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
260      !
261   END FUNCTION sbc_cpl_alloc
262
263
264   SUBROUTINE sbc_cpl_init( k_ice )     
265      !!----------------------------------------------------------------------
266      !!             ***  ROUTINE sbc_cpl_init  ***
267      !!
268      !! ** Purpose :   Initialisation of send and received information from
269      !!                the atmospheric component
270      !!
271      !! ** Method  : * Read namsbc_cpl namelist
272      !!              * define the receive interface
273      !!              * define the send    interface
274      !!              * initialise the OASIS coupler
275      !!----------------------------------------------------------------------
276      INTEGER, INTENT(in) ::   k_ice   ! ice management in the sbc (=0/1/2/3)
277      !
278      INTEGER ::   jn          ! dummy loop index
279      INTEGER ::   ios, inum   ! Local integer
280      REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos
281      !!
282      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             &
283         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  & 
284         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  & 
285         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  & 
286         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  & 
287         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  &
288         &                  sn_rcv_qtr   ,                                                             &
289         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  &
290         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                &
291         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 &
292         &                  sn_rcv_ts_ice, sn_rcv_grnm  , sn_rcv_antm  ,                               &
293         &                  nn_coupled_iceshelf_fluxes  , ln_iceshelf_init_atmos ,                     &
294         &                  rn_greenland_total_fw_flux  , rn_greenland_calving_fraction  ,             &
295         &                  rn_antarctica_total_fw_flux , rn_antarctica_calving_fraction ,             &
296         &                  rn_iceshelf_fluxes_tolerance
297
298                                                                !  Add MEDUSA related fields to namelist
299      NAMELIST/namsbc_cpl/  sn_snd_bio_co2 , sn_snd_bio_dms, sn_snd_bio_chloro,                        &
300         &                  sn_rcv_atm_pco2, sn_rcv_atm_dust
301
302      !!---------------------------------------------------------------------
303      !
304      ! ================================ !
305      !      Namelist informations       !
306      ! ================================ !
307      !
308      IF (ln_timing) CALL timing_start('sbc_cpl_init')
309      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
310      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
311901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist' )
312      !
313      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
314      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
315902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist' )
316      IF(lwm) WRITE ( numond, namsbc_cpl )
317      !
318      IF(lwp) THEN                        ! control print
319         WRITE(numout,*)
320         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
321         WRITE(numout,*)'~~~~~~~~~~~~'
322      ENDIF
323      IF( lwp .AND. ln_cpl ) THEN                        ! control print
324         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
325         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
326         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux
327         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl
328         WRITE(numout,*)'  received fields (mutiple ice categogies)'
329         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
330         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
331         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
332         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
333         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
334         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
335         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
336         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
337         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
338         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
339         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
340         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
341         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')' 
342         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')' 
343         WRITE(numout,*)'      iceberg                         = ', TRIM(sn_rcv_icb%cldes   ), ' (', TRIM(sn_rcv_icb%clcat   ), ')'
344         WRITE(numout,*)'      ice shelf                       = ', TRIM(sn_rcv_isf%cldes   ), ' (', TRIM(sn_rcv_isf%clcat   ), ')'
345         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
346         WRITE(numout,*)'      transmitted solar               = ', TRIM(sn_rcv_qtr%cldes   ), ' (', TRIM(sn_rcv_qtr%clcat   ), ')'
347         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
348         WRITE(numout,*)'      significant wave heigth         = ', TRIM(sn_rcv_hsig%cldes  ), ' (', TRIM(sn_rcv_hsig%clcat  ), ')' 
349         WRITE(numout,*)'      wave to oce energy flux         = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' 
350         WRITE(numout,*)'      Surface Stokes drift grid u     = ', TRIM(sn_rcv_sdrfx%cldes ), ' (', TRIM(sn_rcv_sdrfx%clcat ), ')' 
351         WRITE(numout,*)'      Surface Stokes drift grid v     = ', TRIM(sn_rcv_sdrfy%cldes ), ' (', TRIM(sn_rcv_sdrfy%clcat ), ')' 
352         WRITE(numout,*)'      Mean wave period                = ', TRIM(sn_rcv_wper%cldes  ), ' (', TRIM(sn_rcv_wper%clcat  ), ')' 
353         WRITE(numout,*)'      Mean wave number                = ', TRIM(sn_rcv_wnum%cldes  ), ' (', TRIM(sn_rcv_wnum%clcat  ), ')' 
354         WRITE(numout,*)'      Wave peak frequency             = ', TRIM(sn_rcv_wfreq%cldes ), ' (', TRIM(sn_rcv_wfreq%clcat ), ')'
355         WRITE(numout,*)'      Stress frac adsorbed by waves   = ', TRIM(sn_rcv_tauwoc%cldes), ' (', TRIM(sn_rcv_tauwoc%clcat ), ')' 
356         WRITE(numout,*)'      Stress components by waves      = ', TRIM(sn_rcv_tauw%cldes  ), ' (', TRIM(sn_rcv_tauw%clcat  ), ')'
357         WRITE(numout,*)'      Neutral surf drag coefficient   = ', TRIM(sn_rcv_wdrag%cldes ), ' (', TRIM(sn_rcv_wdrag%clcat ), ')' 
358         WRITE(numout,*)'      Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 
359         WRITE(numout,*)'  sent fields (multiple ice categories)'
360         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
361         WRITE(numout,*)'      top ice layer temperature       = ', TRIM(sn_snd_ttilyr%cldes), ' (', TRIM(sn_snd_ttilyr%clcat), ')'
362         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
363         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
364         WRITE(numout,*)'      total ice fraction              = ', TRIM(sn_snd_ifrac%cldes ), ' (', TRIM(sn_snd_ifrac%clcat ), ')' 
365         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
366         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
367         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
368         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
369         WRITE(numout,*)'      bio co2 flux                    = ', TRIM(sn_snd_bio_co2%cldes),' (', TRIM(sn_snd_bio_co2%clcat), ')'
370         WRITE(numout,*)'      bio dms flux                    = ', TRIM(sn_snd_bio_dms%cldes),' (', TRIM(sn_snd_bio_dms%clcat), ')'
371         WRITE(numout,*)'      bio dms chlorophyll             = ', TRIM(sn_snd_bio_chloro%cldes), ' (', TRIM(sn_snd_bio_chloro%clcat), ')'
372         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
373         WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes  ), ' (', TRIM(sn_snd_cond%clcat  ), ')'
374         WRITE(numout,*)'      meltponds fraction and depth    = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat  ), ')'
375         WRITE(numout,*)'      sea surface freezing temp       = ', TRIM(sn_snd_sstfrz%cldes), ' (', TRIM(sn_snd_sstfrz%clcat), ')'
376         WRITE(numout,*)'      water level                     = ', TRIM(sn_snd_wlev%cldes  ), ' (', TRIM(sn_snd_wlev%clcat  ), ')' 
377         WRITE(numout,*)'      mean sea level pressure         = ', TRIM(sn_rcv_mslp%cldes  ), ' (', TRIM(sn_rcv_mslp%clcat  ), ')' 
378         WRITE(numout,*)'      surface current to waves        = ', TRIM(sn_snd_crtw%cldes  ), ' (', TRIM(sn_snd_crtw%clcat  ), ')' 
379         WRITE(numout,*)'      atm pco2                        = ', TRIM(sn_rcv_atm_pco2%cldes),'(', TRIM(sn_rcv_atm_pco2%clcat   ), ')'
380         WRITE(numout,*)'      atm dust                        = ', TRIM(sn_rcv_atm_dust%cldes),'(', TRIM(sn_rcv_atm_dust%clcat),')'
381         WRITE(numout,*)'                      - referential   = ', sn_snd_crtw%clvref 
382         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor 
383         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd 
384         WRITE(numout,*)'  nn_coupled_iceshelf_fluxes          = ', nn_coupled_iceshelf_fluxes
385         WRITE(numout,*)'  ln_iceshelf_init_atmos              = ', ln_iceshelf_init_atmos
386         WRITE(numout,*)'  rn_greenland_total_fw_flux          = ', rn_greenland_total_fw_flux
387         WRITE(numout,*)'  rn_antarctica_total_fw_flux         = ', rn_antarctica_total_fw_flux
388         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction
389         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction
390         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance
391      ENDIF
392
393      !                                   ! allocate sbccpl arrays
394      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
395     
396      ! ================================ !
397      !   Define the receive interface   !
398      ! ================================ !
399      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
400
401      ! for each field: define the OASIS name                              (srcv(:)%clname)
402      !                 define receive or not from the namelist parameters (srcv(:)%laction)
403      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
404
405      ! default definitions of srcv
406      srcv(:)%laction = .FALSE. 
407      srcv(:)%clgrid = 'T' 
408      srcv(:)%nsgn = 1. 
409      srcv(:)%nct = 1 
410      srcv(:)%dimensions = 2 
411
412      !                                                      ! ------------------------- !
413      !                                                      ! ice and ocean wind stress !   
414      !                                                      ! ------------------------- !
415      !                                                           ! Name
416      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
417      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
418      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
419      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
420      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
421      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
422      !
423      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
424      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
425      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
426      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
427      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
428      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
429      !
430      ! Vectors: change of sign at north fold ONLY if on the local grid
431      IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  &
432           .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled
433      !
434      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
435     
436      !                                                           ! Set grid and action
437      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
438      CASE( 'T' ) 
439         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
440         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
441         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
442      CASE( 'U,V' ) 
443         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
444         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
445         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
446         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
447         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
448      CASE( 'U,V,T' )
449         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
450         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
451         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
452         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
453         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
454      CASE( 'U,V,I' )
455         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
456         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
457         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
458         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
459         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
460      CASE( 'U,V,F' )
461         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
462         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
463         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
464         !srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
465! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment
466         srcv(jpr_otx1)%laction = .TRUE. 
467         srcv(jpr_oty1)%laction = .TRUE.
468!
469         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
470      CASE( 'T,I' ) 
471         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
472         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
473         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
474         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
475      CASE( 'T,F' ) 
476         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
477         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
478         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
479         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
480      CASE( 'T,U,V' )
481         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
482         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
483         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
484         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
485         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
486      CASE default   
487         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
488      END SELECT
489      !
490      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
491         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
492      !
493      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
494            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
495            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
496            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
497            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
498      ENDIF
499      !
500      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
501         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
502         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
503         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
504      ENDIF
505      ENDIF
506
507      !                                                      ! ------------------------- !
508      !                                                      !    freshwater budget      !   E-P
509      !                                                      ! ------------------------- !
510      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
511      ! over ice of free ocean within the same atmospheric cell.cd
512      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
513      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
514      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
515      srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation
516      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
517      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
518      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
519      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
520      CASE( 'none'          )       ! nothing to do
521      CASE( 'oce only'      )   ;   srcv(jpr_oemp)%laction = .TRUE. 
522      CASE( 'conservative'  )
523         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
524         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
525      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
526      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
527      END SELECT
528      !
529      !                                                      ! ------------------------- !
530      !                                                      !     Runoffs & Calving     !   
531      !                                                      ! ------------------------- !
532      srcv(jpr_rnf   )%clname = 'O_Runoff'
533      srcv(jpr_rnf_1d   )%clname = 'ORunff1D' 
534      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' .OR. TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN 
535         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) srcv(jpr_rnf)%laction = .TRUE. 
536         IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) THEN
537            srcv(jpr_rnf_1d)%laction = .TRUE. 
538            srcv(jpr_rnf_1d)%dimensions = 1 ! 1D field passed through coupler
539         END IF
540         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
541         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
542         IF(lwp) WRITE(numout,*)
543         IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
544      ENDIF
545      !
546      srcv(jpr_cal   )%clname = 'OCalving'   
547      IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.     
548 
549      srcv(jpr_grnm  )%clname = 'OGrnmass' 
550      IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' ) THEN
551         srcv(jpr_grnm)%laction = .TRUE. 
552         srcv(jpr_grnm)%dimensions = 0 ! Scalar field
553      ENDIF
554       
555      srcv(jpr_antm  )%clname = 'OAntmass' 
556      IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' ) THEN
557         srcv(jpr_antm)%laction = .TRUE.
558         srcv(jpr_antm)%dimensions = 0 ! Scalar field
559      ENDIF
560
561      srcv(jpr_isf)%clname = 'OIcshelf'   ;  IF( TRIM( sn_rcv_isf%cldes) == 'coupled' )   srcv(jpr_isf)%laction = .TRUE.
562      srcv(jpr_icb)%clname = 'OIceberg'   ;  IF( TRIM( sn_rcv_icb%cldes) == 'coupled' )   srcv(jpr_icb)%laction = .TRUE.
563
564      IF( srcv(jpr_isf)%laction .AND. ln_isf ) THEN
565         l_isfcpl             = .TRUE.                      ! -> no need to read isf in sbcisf
566         IF(lwp) WRITE(numout,*)
567         IF(lwp) WRITE(numout,*) '   iceshelf received from oasis '
568      ENDIF
569      !
570      !                                                      ! ------------------------- !
571      !                                                      !    non solar radiation    !   Qns
572      !                                                      ! ------------------------- !
573      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
574      srcv(jpr_qnsice)%clname = 'O_QnsIce'
575      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
576      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
577      CASE( 'none'          )       ! nothing to do
578      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
579      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
580      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
581      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
582      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
583      END SELECT
584      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. nn_cats_cpl > 1 ) &
585         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
586      !
587      !                                                      ! ------------------------- !
588      !                                                      !    solar radiation        !   Qsr
589      !                                                      ! ------------------------- !
590      srcv(jpr_qsroce)%clname = 'O_QsrOce'
591      srcv(jpr_qsrice)%clname = 'O_QsrIce'
592      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
593      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
594      CASE( 'none'          )       ! nothing to do
595      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
596      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
597      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
598      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
599      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
600      END SELECT
601      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. nn_cats_cpl > 1 ) &
602         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
603      !
604      !                                                      ! ------------------------- !
605      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
606      !                                                      ! ------------------------- !
607      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
608      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
609      !
610      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
611      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' )  &
612         &   CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
613      !
614      !                                                      ! ------------------------- !
615      !                                                      !      10m wind module      !   
616      !                                                      ! ------------------------- !
617      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
618      !
619      !                                                      ! ------------------------- !
620      !                                                      !   wind stress module      !   
621      !                                                      ! ------------------------- !
622      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
623      lhftau = srcv(jpr_taum)%laction
624      !
625      !                                                      ! ------------------------- !
626      !                                                      !      Atmospheric CO2      !
627      !                                                      ! ------------------------- !
628      srcv(jpr_co2 )%clname = 'O_AtmCO2'   
629      IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )  THEN
630         srcv(jpr_co2 )%laction = .TRUE.
631         l_co2cpl = .TRUE.
632         IF(lwp) WRITE(numout,*)
633         IF(lwp) WRITE(numout,*) '   Atmospheric pco2 received from oasis '
634         IF(lwp) WRITE(numout,*)
635      ENDIF
636      !                                                      ! --------------------------------------- !
637      !                                                      ! Incoming CO2 and DUST fluxes for MEDUSA !
638      !                                                      ! --------------------------------------- !
639      srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
640      IF (TRIM(sn_rcv_atm_pco2%cldes) == 'medusa') THEN
641        srcv(jpr_atm_pco2)%laction = .TRUE.
642      END IF
643
644      srcv(jpr_atm_dust)%clname = 'OATMDUST'
645      IF (TRIM(sn_rcv_atm_dust%cldes) == 'medusa')  THEN
646        srcv(jpr_atm_dust)%laction = .TRUE.
647      END IF
648      !
649      !                                                      ! ------------------------- !
650      !                                                      ! Mean Sea Level Pressure   !
651      !                                                      ! ------------------------- !
652      srcv(jpr_mslp)%clname = 'O_MSLP'     ;   IF( TRIM(sn_rcv_mslp%cldes  ) == 'coupled' )    srcv(jpr_mslp)%laction = .TRUE. 
653      !
654      !                                                      ! ------------------------- !
655      !                                                      !  ice topmelt and botmelt  !   
656      !                                                      ! ------------------------- !
657      srcv(jpr_topm )%clname = 'OTopMlt'
658      srcv(jpr_botm )%clname = 'OBotMlt'
659      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
660         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
661            srcv(jpr_topm:jpr_botm)%nct = nn_cats_cpl
662         ELSE
663            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
664         ENDIF
665         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
666      ENDIF
667      !                                                      ! ------------------------- !
668      !                                                      !    transmitted solar      !   
669      !                                                      ! ------------------------- !
670      srcv(jpr_qtr )%clname = 'OQtr'
671      IF( TRIM(sn_rcv_qtr%cldes) == 'coupled' ) THEN
672         IF ( TRIM( sn_rcv_qtr%clcat ) == 'yes' ) THEN
673            srcv(jpr_qtr)%nct = nn_cats_cpl
674         ELSE
675            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qtr%clcat should always be set to yes currently' )
676         ENDIF
677         srcv(jpr_qtr)%laction = .TRUE.
678      ENDIF
679
680      !                                                      ! ------------------------- !
681      !                                                      !    ice skin temperature   !   
682      !                                                      ! ------------------------- !
683      srcv(jpr_ts_ice)%clname = 'OTsfIce'    ! needed by Met Office
684      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' )   srcv(jpr_ts_ice)%laction = .TRUE.
685      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' )   srcv(jpr_ts_ice)%nct     = nn_cats_cpl
686      IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl
687
688#if defined key_si3
689      IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
690         IF( .NOT.srcv(jpr_ts_ice)%laction )  &
691            &   CALL ctl_stop( 'sbc_cpl_init: srcv(jpr_ts_ice)%laction should be set to true when ln_cndflx=T' )     
692      ENDIF
693#endif
694      !                                                      ! ------------------------- !
695      !                                                      !      Wave breaking        !   
696      !                                                      ! ------------------------- !
697      srcv(jpr_hsig)%clname  = 'O_Hsigwa'    ! significant wave height
698      IF( TRIM(sn_rcv_hsig%cldes  ) == 'coupled' )  THEN
699         srcv(jpr_hsig)%laction = .TRUE.
700         cpl_hsig = .TRUE.
701      ENDIF
702      srcv(jpr_phioc)%clname = 'O_PhiOce'    ! wave to ocean energy
703      IF( TRIM(sn_rcv_phioc%cldes ) == 'coupled' )  THEN
704         srcv(jpr_phioc)%laction = .TRUE.
705         cpl_phioc = .TRUE.
706      ENDIF
707      srcv(jpr_sdrftx)%clname = 'O_Sdrfx'    ! Stokes drift in the u direction
708      IF( TRIM(sn_rcv_sdrfx%cldes ) == 'coupled' )  THEN
709         srcv(jpr_sdrftx)%laction = .TRUE.
710         cpl_sdrftx = .TRUE.
711      ENDIF
712      srcv(jpr_sdrfty)%clname = 'O_Sdrfy'    ! Stokes drift in the v direction
713      IF( TRIM(sn_rcv_sdrfy%cldes ) == 'coupled' )  THEN
714         srcv(jpr_sdrfty)%laction = .TRUE.
715         cpl_sdrfty = .TRUE.
716      ENDIF
717      srcv(jpr_wper)%clname = 'O_WPer'       ! mean wave period
718      IF( TRIM(sn_rcv_wper%cldes  ) == 'coupled' )  THEN
719         srcv(jpr_wper)%laction = .TRUE.
720         cpl_wper = .TRUE.
721      ENDIF
722      srcv(jpr_wfreq)%clname = 'O_WFreq'     ! wave peak frequency
723      IF( TRIM(sn_rcv_wfreq%cldes ) == 'coupled' )  THEN
724         srcv(jpr_wfreq)%laction = .TRUE.
725         cpl_wfreq = .TRUE.
726      ENDIF
727      srcv(jpr_wnum)%clname = 'O_WNum'       ! mean wave number
728      IF( TRIM(sn_rcv_wnum%cldes ) == 'coupled' )  THEN
729         srcv(jpr_wnum)%laction = .TRUE.
730         cpl_wnum = .TRUE.
731      ENDIF
732      srcv(jpr_tauwoc)%clname = 'O_TauOce'   ! stress fraction adsorbed by the wave
733      IF( TRIM(sn_rcv_tauwoc%cldes ) == 'coupled' )  THEN
734         srcv(jpr_tauwoc)%laction = .TRUE.
735         cpl_tauwoc = .TRUE.
736      ENDIF
737      srcv(jpr_tauwx)%clname = 'O_Tauwx'      ! ocean stress from wave in the x direction
738      srcv(jpr_tauwy)%clname = 'O_Tauwy'      ! ocean stress from wave in the y direction
739      IF( TRIM(sn_rcv_tauw%cldes ) == 'coupled' )  THEN
740         srcv(jpr_tauwx)%laction = .TRUE.
741         srcv(jpr_tauwy)%laction = .TRUE.
742         cpl_tauw = .TRUE.
743      ENDIF
744      srcv(jpr_wdrag)%clname = 'O_WDrag'     ! neutral surface drag coefficient
745      IF( TRIM(sn_rcv_wdrag%cldes ) == 'coupled' )  THEN
746         srcv(jpr_wdrag)%laction = .TRUE.
747         cpl_wdrag = .TRUE.
748      ENDIF
749      IF( srcv(jpr_tauwoc)%laction .AND. srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction ) &
750            CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
751                                     '(sn_rcv_tauwoc=coupled and sn_rcv_tauw=coupled)' )
752      !
753      !                                                      ! ------------------------------- !
754      !                                                      !   OPA-SAS coupling - rcv by opa !   
755      !                                                      ! ------------------------------- !
756      srcv(jpr_sflx)%clname = 'O_SFLX'
757      srcv(jpr_fice)%clname = 'RIceFrc'
758      !
759      IF( nn_components == jp_iam_opa ) THEN    ! OPA coupled to SAS via OASIS: force received field by OPA (sent by SAS)
760         srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
761         srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
762         srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
763         srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
764         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
765         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
766         ! Vectors: change of sign at north fold ONLY if on the local grid
767         srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
768         sn_rcv_tau%clvgrd = 'U,V'
769         sn_rcv_tau%clvor = 'local grid'
770         sn_rcv_tau%clvref = 'spherical'
771         sn_rcv_emp%cldes = 'oce only'
772         !
773         IF(lwp) THEN                        ! control print
774            WRITE(numout,*)
775            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
776            WRITE(numout,*)'               OPA component  '
777            WRITE(numout,*)
778            WRITE(numout,*)'  received fields from SAS component '
779            WRITE(numout,*)'                  ice cover '
780            WRITE(numout,*)'                  oce only EMP  '
781            WRITE(numout,*)'                  salt flux  '
782            WRITE(numout,*)'                  mixed oce-ice solar flux  '
783            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
784            WRITE(numout,*)'                  wind stress U,V on local grid and sperical coordinates '
785            WRITE(numout,*)'                  wind stress module'
786            WRITE(numout,*)
787         ENDIF
788      ENDIF
789      !                                                      ! -------------------------------- !
790      !                                                      !   OPA-SAS coupling - rcv by sas  !   
791      !                                                      ! -------------------------------- !
792      srcv(jpr_toce  )%clname = 'I_SSTSST'
793      srcv(jpr_soce  )%clname = 'I_SSSal'
794      srcv(jpr_ocx1  )%clname = 'I_OCurx1'
795      srcv(jpr_ocy1  )%clname = 'I_OCury1'
796      srcv(jpr_ssh   )%clname = 'I_SSHght'
797      srcv(jpr_e3t1st)%clname = 'I_E3T1st'   
798      srcv(jpr_fraqsr)%clname = 'I_FraQsr'   
799      !
800      IF( nn_components == jp_iam_sas ) THEN
801         IF( .NOT. ln_cpl ) srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
802         IF( .NOT. ln_cpl ) srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
803         IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
804         srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
805         srcv( jpr_e3t1st )%laction = .NOT.ln_linssh
806         srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
807         srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
808         ! Vectors: change of sign at north fold ONLY if on the local grid
809         srcv(jpr_ocx1:jpr_ocy1)%nsgn = -1.
810         ! Change first letter to couple with atmosphere if already coupled OPA
811         ! this is nedeed as each variable name used in the namcouple must be unique:
812         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
813         DO jn = 1, jprcv
814            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
815         END DO
816         !
817         IF(lwp) THEN                        ! control print
818            WRITE(numout,*)
819            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
820            WRITE(numout,*)'               SAS component  '
821            WRITE(numout,*)
822            IF( .NOT. ln_cpl ) THEN
823               WRITE(numout,*)'  received fields from OPA component '
824            ELSE
825               WRITE(numout,*)'  Additional received fields from OPA component : '
826            ENDIF
827            WRITE(numout,*)'               sea surface temperature (Celsius) '
828            WRITE(numout,*)'               sea surface salinity ' 
829            WRITE(numout,*)'               surface currents ' 
830            WRITE(numout,*)'               sea surface height ' 
831            WRITE(numout,*)'               thickness of first ocean T level '       
832            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
833            WRITE(numout,*)
834         ENDIF
835      ENDIF
836
837      ! ================================ !
838      !     Define the send interface    !
839      ! ================================ !
840      ! for each field: define the OASIS name                           (ssnd(:)%clname)
841      !                 define send or not from the namelist parameters (ssnd(:)%laction)
842      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
843     
844      ! default definitions of nsnd
845      ssnd(:)%laction = .FALSE. 
846      ssnd(:)%clgrid = 'T' 
847      ssnd(:)%nsgn = 1. 
848      ssnd(:)%nct = 1 
849      ssnd(:)%dimensions = 2 
850         
851      !                                                      ! ------------------------- !
852      !                                                      !    Surface temperature    !
853      !                                                      ! ------------------------- !
854      ssnd(jps_toce)%clname   = 'O_SSTSST'
855      ssnd(jps_tice)%clname   = 'O_TepIce'
856      ssnd(jps_ttilyr)%clname = 'O_TtiLyr'
857      ssnd(jps_tmix)%clname   = 'O_TepMix'
858      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
859      CASE( 'none'                                 )       ! nothing to do
860      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
861      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice' )
862         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
863         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = nn_cats_cpl
864      CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE.
865      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
866      END SELECT
867           
868      !                                                      ! ------------------------- !
869      !                                                      !          Albedo           !
870      !                                                      ! ------------------------- !
871      ssnd(jps_albice)%clname = 'O_AlbIce' 
872      ssnd(jps_albmix)%clname = 'O_AlbMix'
873      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
874      CASE( 'none'                 )     ! nothing to do
875      CASE( 'ice' , 'weighted ice' )   ; ssnd(jps_albice)%laction = .TRUE.
876      CASE( 'mixed oce-ice'        )   ; ssnd(jps_albmix)%laction = .TRUE.
877      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
878      END SELECT
879      !
880      ! Need to calculate oceanic albedo if
881      !     1. sending mixed oce-ice albedo or
882      !     2. receiving mixed oce-ice solar radiation
883      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
884         CALL oce_alb( zaos, zacs )
885         ! Due to lack of information on nebulosity : mean clear/overcast sky
886         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
887      ENDIF
888      !                                                      ! ------------------------- !
889      !                                                      !  Ice fraction & Thickness !
890      !                                                      ! ------------------------- !
891      ssnd(jps_fice)%clname  = 'OIceFrc'
892      ssnd(jps_ficet)%clname = 'OIceFrcT' 
893      ssnd(jps_hice)%clname  = 'OIceTck'
894      ssnd(jps_a_p)%clname   = 'OPndFrc'
895      ssnd(jps_ht_p)%clname  = 'OPndTck'
896      ssnd(jps_hsnw)%clname  = 'OSnwTck'
897      ssnd(jps_fice1)%clname = 'OIceFrd'
898      IF( k_ice /= 0 ) THEN
899         ssnd(jps_fice)%laction  = .TRUE.                 ! if ice treated in the ocean (even in climato case)
900         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used producing atmos-to-ice fluxes (Met Office requirement)
901! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
902         IF ( TRIM( sn_snd_thick%clcat  ) == 'yes' ) ssnd(jps_fice)%nct  = nn_cats_cpl
903         IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = nn_cats_cpl
904      ENDIF
905     
906      IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE. 
907
908      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
909      CASE( 'none'         )       ! nothing to do
910      CASE( 'ice and snow' ) 
911         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
912         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
913            ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl
914         ENDIF
915      CASE ( 'weighted ice and snow' ) 
916         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
917         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl
918      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
919      END SELECT
920
921      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office)
922#if defined key_si3 || defined key_cice
923       a_i_last_couple(:,:,:) = 0._wp
924#endif
925      !                                                      ! ------------------------- !
926      !                                                      !      Ice Meltponds        !
927      !                                                      ! ------------------------- !
928      ! Needed by Met Office
929      ssnd(jps_a_p)%clname  = 'OPndFrc'   
930      ssnd(jps_ht_p)%clname = 'OPndTck'   
931      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) ) 
932      CASE ( 'none' ) 
933         ssnd(jps_a_p)%laction  = .FALSE. 
934         ssnd(jps_ht_p)%laction = .FALSE. 
935      CASE ( 'ice only' ) 
936         ssnd(jps_a_p)%laction  = .TRUE. 
937         ssnd(jps_ht_p)%laction = .TRUE. 
938         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
939            ssnd(jps_a_p)%nct  = nn_cats_cpl 
940            ssnd(jps_ht_p)%nct = nn_cats_cpl 
941         ELSE
942            IF ( nn_cats_cpl > 1 ) THEN
943               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' ) 
944            ENDIF
945         ENDIF
946      CASE ( 'weighted ice' ) 
947         ssnd(jps_a_p)%laction  = .TRUE. 
948         ssnd(jps_ht_p)%laction = .TRUE. 
949         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
950            ssnd(jps_a_p)%nct  = nn_cats_cpl 
951            ssnd(jps_ht_p)%nct = nn_cats_cpl 
952         ENDIF
953      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes; '//sn_snd_mpnd%cldes ) 
954      END SELECT 
955 
956      !                                                      ! ------------------------- !
957      !                                                      !      Surface current      !
958      !                                                      ! ------------------------- !
959      !        ocean currents              !            ice velocities
960      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
961      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
962      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
963      ssnd(jps_ocxw)%clname = 'O_OCurxw' 
964      ssnd(jps_ocyw)%clname = 'O_OCuryw' 
965      !
966      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
967
968      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
969         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
970      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
971         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
972      ENDIF
973      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
974      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
975      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
976      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
977      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
978      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
979      CASE( 'weighted oce and ice' )   !   nothing to do
980      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
981      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
982      END SELECT
983
984      ssnd(jps_ocxw:jps_ocyw)%nsgn = -1.   ! vectors: change of the sign at the north fold
985       
986      IF( sn_snd_crtw%clvgrd == 'U,V' ) THEN
987         ssnd(jps_ocxw)%clgrid = 'U' ; ssnd(jps_ocyw)%clgrid = 'V' 
988      ELSE IF( sn_snd_crtw%clvgrd /= 'T' ) THEN
989         CALL ctl_stop( 'sn_snd_crtw%clvgrd must be equal to T' ) 
990      ENDIF
991      IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) ssnd(jps_ocxw:jps_ocyw)%nsgn = 1. 
992      SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
993         CASE( 'none'                 )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .FALSE. 
994         CASE( 'oce only'             )   ; ssnd(jps_ocxw:jps_ocyw)%laction = .TRUE. 
995         CASE( 'weighted oce and ice' )   !   nothing to do
996         CASE( 'mixed oce-ice'        )   ; ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE. 
997         CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crtw%cldes' ) 
998      END SELECT 
999
1000      !                                                      ! ------------------------- !
1001      !                                                      !          CO2 flux         !
1002      !                                                      ! ------------------------- !
1003      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
1004      !
1005      !                                                      ! ------------------------- !
1006      !                                                      !   MEDUSA output fields    !
1007      !                                                      ! ------------------------- !
1008      ! Surface dimethyl sulphide from Medusa
1009      ssnd(jps_bio_dms)%clname = 'OBioDMS'
1010      IF( TRIM(sn_snd_bio_dms%cldes) == 'medusa' )    ssnd(jps_bio_dms )%laction = .TRUE.
1011
1012      ! Surface CO2 flux from Medusa
1013      ssnd(jps_bio_co2)%clname = 'OBioCO2'
1014      IF( TRIM(sn_snd_bio_co2%cldes) == 'medusa' )    ssnd(jps_bio_co2 )%laction = .TRUE.
1015
1016      ! Surface chlorophyll from Medusa
1017      ssnd(jps_bio_chloro)%clname = 'OBioChlo'
1018      IF( TRIM(sn_snd_bio_chloro%cldes) == 'medusa' ) ssnd(jps_bio_chloro )%laction = .TRUE.
1019      !
1020      !                                                      ! ------------------------- !
1021      !                                                      ! Sea surface freezing temp !
1022      !                                                      ! ------------------------- !
1023      ! needed by Met Office
1024      ssnd(jps_sstfrz)%clname = 'O_SSTFrz' ; IF( TRIM(sn_snd_sstfrz%cldes) == 'coupled' )  ssnd(jps_sstfrz)%laction = .TRUE. 
1025      !
1026      !                                                      ! ------------------------- !
1027      !                                                      !    Ice conductivity       !
1028      !                                                      ! ------------------------- !
1029      ! needed by Met Office
1030      ! Note that ultimately we will move to passing an ocean effective conductivity as well so there
1031      ! will be some changes to the parts of the code which currently relate only to ice conductivity
1032      ssnd(jps_ttilyr )%clname = 'O_TtiLyr' 
1033      SELECT CASE ( TRIM( sn_snd_ttilyr%cldes ) ) 
1034      CASE ( 'none' ) 
1035         ssnd(jps_ttilyr)%laction = .FALSE. 
1036      CASE ( 'ice only' ) 
1037         ssnd(jps_ttilyr)%laction = .TRUE. 
1038         IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN
1039            ssnd(jps_ttilyr)%nct = nn_cats_cpl 
1040         ELSE
1041            IF ( nn_cats_cpl > 1 ) THEN
1042               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_ttilyr%cldes if not exchanging category fields' ) 
1043            ENDIF
1044         ENDIF
1045      CASE ( 'weighted ice' ) 
1046         ssnd(jps_ttilyr)%laction = .TRUE. 
1047         IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl 
1048      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_ttilyr%cldes;'//sn_snd_ttilyr%cldes ) 
1049      END SELECT
1050
1051      ssnd(jps_kice )%clname = 'OIceKn' 
1052      SELECT CASE ( TRIM( sn_snd_cond%cldes ) ) 
1053      CASE ( 'none' ) 
1054         ssnd(jps_kice)%laction = .FALSE. 
1055      CASE ( 'ice only' ) 
1056         ssnd(jps_kice)%laction = .TRUE. 
1057         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN
1058            ssnd(jps_kice)%nct = nn_cats_cpl 
1059         ELSE
1060            IF ( nn_cats_cpl > 1 ) THEN
1061               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' ) 
1062            ENDIF
1063         ENDIF
1064      CASE ( 'weighted ice' ) 
1065         ssnd(jps_kice)%laction = .TRUE. 
1066         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl 
1067      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes;'//sn_snd_cond%cldes ) 
1068      END SELECT 
1069      !
1070      !                                                      ! ------------------------- !
1071      !                                                      !     Sea surface height    !
1072      !                                                      ! ------------------------- !
1073      ssnd(jps_wlev)%clname = 'O_Wlevel' ;  IF( TRIM(sn_snd_wlev%cldes) == 'coupled' )   ssnd(jps_wlev)%laction = .TRUE. 
1074
1075      !                                                      ! ------------------------------- !
1076      !                                                      !   OPA-SAS coupling - snd by opa !   
1077      !                                                      ! ------------------------------- !
1078      ssnd(jps_ssh   )%clname = 'O_SSHght' 
1079      ssnd(jps_soce  )%clname = 'O_SSSal' 
1080      ssnd(jps_e3t1st)%clname = 'O_E3T1st'   
1081      ssnd(jps_fraqsr)%clname = 'O_FraQsr'
1082      !
1083      IF( nn_components == jp_iam_opa ) THEN
1084         ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
1085         ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
1086         ssnd( jps_e3t1st )%laction = .NOT.ln_linssh
1087         ! vector definition: not used but cleaner...
1088         ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
1089         ssnd(jps_ocy1)%clgrid  = 'V'        !           and           V-point
1090         sn_snd_crt%clvgrd = 'U,V'
1091         sn_snd_crt%clvor = 'local grid'
1092         sn_snd_crt%clvref = 'spherical'
1093         !
1094         IF(lwp) THEN                        ! control print
1095            WRITE(numout,*)
1096            WRITE(numout,*)'  sent fields to SAS component '
1097            WRITE(numout,*)'               sea surface temperature (T before, Celsius) '
1098            WRITE(numout,*)'               sea surface salinity ' 
1099            WRITE(numout,*)'               surface currents U,V on local grid and spherical coordinates' 
1100            WRITE(numout,*)'               sea surface height ' 
1101            WRITE(numout,*)'               thickness of first ocean T level '       
1102            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
1103            WRITE(numout,*)
1104         ENDIF
1105      ENDIF
1106      !                                                      ! ------------------------------- !
1107      !                                                      !   OPA-SAS coupling - snd by sas !   
1108      !                                                      ! ------------------------------- !
1109      ssnd(jps_sflx  )%clname = 'I_SFLX'     
1110      ssnd(jps_fice2 )%clname = 'IIceFrc'
1111      ssnd(jps_qsroce)%clname = 'I_QsrOce'   
1112      ssnd(jps_qnsoce)%clname = 'I_QnsOce'   
1113      ssnd(jps_oemp  )%clname = 'IOEvaMPr' 
1114      ssnd(jps_otx1  )%clname = 'I_OTaux1'   
1115      ssnd(jps_oty1  )%clname = 'I_OTauy1'   
1116      ssnd(jps_rnf   )%clname = 'I_Runoff'   
1117      ssnd(jps_taum  )%clname = 'I_TauMod'   
1118      !
1119      IF( nn_components == jp_iam_sas ) THEN
1120         IF( .NOT. ln_cpl ) ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
1121         ssnd( (/jps_qsroce, jps_qnsoce, jps_oemp, jps_fice2, jps_sflx, jps_otx1, jps_oty1, jps_taum/) )%laction = .TRUE.
1122         !
1123         ! Change first letter to couple with atmosphere if already coupled with sea_ice
1124         ! this is nedeed as each variable name used in the namcouple must be unique:
1125         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
1126         DO jn = 1, jpsnd
1127            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
1128         END DO
1129         !
1130         IF(lwp) THEN                        ! control print
1131            WRITE(numout,*)
1132            IF( .NOT. ln_cpl ) THEN
1133               WRITE(numout,*)'  sent fields to OPA component '
1134            ELSE
1135               WRITE(numout,*)'  Additional sent fields to OPA component : '
1136            ENDIF
1137            WRITE(numout,*)'                  ice cover '
1138            WRITE(numout,*)'                  oce only EMP  '
1139            WRITE(numout,*)'                  salt flux  '
1140            WRITE(numout,*)'                  mixed oce-ice solar flux  '
1141            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
1142            WRITE(numout,*)'                  wind stress U,V components'
1143            WRITE(numout,*)'                  wind stress module'
1144         ENDIF
1145      ENDIF
1146
1147      ! Initialise 1D river outflow scheme
1148      nn_cpl_river = 1 
1149      IF ( TRIM( sn_rcv_rnf%cldes ) == 'coupled1d' ) CALL cpl_rnf_1d_init   ! Coupled runoff using 1D array
1150     
1151      ! =================================================== !
1152      ! Allocate all parts of frcv used for received fields !
1153      ! =================================================== !
1154      DO jn = 1, jprcv
1155
1156         IF ( srcv(jn)%laction ) THEN
1157            SELECT CASE( srcv(jn)%dimensions )
1158            !
1159            CASE( 0 )   ! Scalar field
1160               ALLOCATE( frcv(jn)%z3(1,1,1) )
1161               
1162            CASE( 1 )   ! 1D field
1163               ALLOCATE( frcv(jn)%z3(nn_cpl_river,1,1) )
1164               
1165            CASE DEFAULT  ! 2D (or pseudo 3D) field.
1166               ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
1167               
1168            END SELECT
1169         END IF
1170
1171      END DO
1172      ! Allocate taum part of frcv which is used even when not received as coupling field
1173      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
1174      ! Allocate w10m part of frcv which is used even when not received as coupling field
1175      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
1176      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
1177      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
1178      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
1179      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
1180      IF( k_ice /= 0 ) THEN
1181         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
1182         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
1183      END IF
1184
1185      !
1186      ! ================================ !
1187      !   initialisation of the coupler  !
1188      ! ================================ !
1189
1190      CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
1191     
1192      IF (ln_usecplmask) THEN
1193         xcplmask(:,:,:) = 0.
1194         CALL iom_open( 'cplmask', inum )
1195         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
1196            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
1197         CALL iom_close( inum )
1198      ELSE
1199         xcplmask(:,:,:) = 1.
1200      ENDIF
1201      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
1202      !
1203      IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN 
1204          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
1205          ! more complicated could be done if required.
1206          greenland_icesheet_mask = 0.0 
1207          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0 
1208          antarctica_icesheet_mask = 0.0 
1209          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0 
1210 
1211          IF( .not. ln_rstart ) THEN
1212             greenland_icesheet_mass = 0.0 
1213             greenland_icesheet_mass_rate_of_change = 0.0 
1214             greenland_icesheet_timelapsed = 0.0 
1215             antarctica_icesheet_mass = 0.0 
1216             antarctica_icesheet_mass_rate_of_change = 0.0 
1217             antarctica_icesheet_timelapsed = 0.0 
1218          ENDIF
1219 
1220      ENDIF 
1221      !
1222      IF (ln_timing) CALL timing_stop('sbc_cpl_init')
1223      !
1224   END SUBROUTINE sbc_cpl_init
1225
1226
1227   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
1228      !!----------------------------------------------------------------------
1229      !!             ***  ROUTINE sbc_cpl_rcv  ***
1230      !!
1231      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
1232      !!                provide the ocean heat and freshwater fluxes.
1233      !!
1234      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
1235      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
1236      !!                to know if the field was really received or not
1237      !!
1238      !!              --> If ocean stress was really received:
1239      !!
1240      !!                  - transform the received ocean stress vector from the received
1241      !!                 referential and grid into an atmosphere-ocean stress in
1242      !!                 the (i,j) ocean referencial and at the ocean velocity point.
1243      !!                    The received stress are :
1244      !!                     - defined by 3 components (if cartesian coordinate)
1245      !!                            or by 2 components (if spherical)
1246      !!                     - oriented along geographical   coordinate (if eastward-northward)
1247      !!                            or  along the local grid coordinate (if local grid)
1248      !!                     - given at U- and V-point, resp.   if received on 2 grids
1249      !!                            or at T-point               if received on 1 grid
1250      !!                    Therefore and if necessary, they are successively
1251      !!                  processed in order to obtain them
1252      !!                     first  as  2 components on the sphere
1253      !!                     second as  2 components oriented along the local grid
1254      !!                     third  as  2 components on the U,V grid
1255      !!
1256      !!              -->
1257      !!
1258      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
1259      !!             and total ocean freshwater fluxes 
1260      !!
1261      !! ** Method  :   receive all fields from the atmosphere and transform
1262      !!              them into ocean surface boundary condition fields
1263      !!
1264      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
1265      !!                        taum         wind stress module at T-point
1266      !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
1267      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
1268      !!                                     and the latent heat flux of solid precip. melting
1269      !!                        qsr          solar ocean heat fluxes   (ocean only case)
1270      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
1271      !!----------------------------------------------------------------------
1272      USE zdf_oce,  ONLY :   ln_zdfswm
1273      !
1274      INTEGER, INTENT(in) ::   kt          ! ocean model time step index
1275      INTEGER, INTENT(in) ::   k_fsbc      ! frequency of sbc (-> ice model) computation
1276      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
1277      !!
1278      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
1279      INTEGER  ::   ji, jj, jn             ! dummy loop indices
1280      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
1281      INTEGER  ::   ikchoix
1282      REAL(wp), DIMENSION(jpi,jpj) ::   ztx2, zty2
1283      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
1284      REAL(wp) ::   zcoef                  ! temporary scalar
1285      LOGICAL  ::   ll_wrtstp              ! write diagnostics?
1286      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
1287      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
1288      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in 
1289      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b 
1290      REAL(wp) ::   zmask_sum, zepsilon   
1291      REAL(wp) ::   zzx, zzy               ! temporary variables
1292      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra
1293      !!----------------------------------------------------------------------
1294      !
1295      !
1296      IF (ln_timing) CALL timing_start('sbc_cpl_rcv')
1297      !
1298      ll_wrtstp  = (( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )) .AND. (nn_print>0)
1299      !
1300      IF( kt == nit000 ) THEN
1301      !   cannot be done in the init phase when we use agrif as cpl_freq requires that oasis_enddef is done
1302         ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
1303         IF( ln_dm2dc .AND. ncpl_qsr_freq /= 86400 )   &
1304            &   CALL ctl_stop( 'sbc_cpl_rcv: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
1305         ncpl_qsr_freq = 86400 / ncpl_qsr_freq   ! used by top
1306      ENDIF
1307      !
1308      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1309      !
1310      !                                                      ! ======================================================= !
1311      !                                                      ! Receive all the atmos. fields (including ice information)
1312      !                                                      ! ======================================================= !
1313      isec = ( kt - nit000 ) * NINT( rdt )                      ! date of exchanges
1314      DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
1315        IF( srcv(jn)%laction ) THEN 
1316 
1317          IF ( srcv(jn)%dimensions <= 1 ) THEN
1318            CALL cpl_rcv_1d( jn, isec, frcv(jn)%z3, SIZE(frcv(jn)%z3), nrcvinfo(jn) ) 
1319          ELSE
1320            CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) ) 
1321          END IF
1322
1323        END IF
1324      END DO
1325
1326      !                                                      ! ========================= !
1327      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
1328         !                                                   ! ========================= !
1329         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
1330         ! => need to be done only when we receive the field
1331         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
1332            !
1333            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1334               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1335               !
1336               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
1337                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
1338               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1339               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1340               !
1341               IF( srcv(jpr_otx2)%laction ) THEN
1342                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
1343                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
1344                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1345                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1346               ENDIF
1347               !
1348            ENDIF
1349            !
1350            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1351               !                                                       ! (geographical to local grid -> rotate the components)
1352               IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN
1353                  ! Temporary code for HadGEM3 - will be removed eventually.
1354        ! Only applies when we have only taux on U grid and tauy on V grid
1355             DO jj=2,jpjm1
1356                DO ji=2,jpim1
1357                     ztx(ji,jj)=0.25*vmask(ji,jj,1)                &
1358                        *(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1)    &
1359                        +frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))
1360                     zty(ji,jj)=0.25*umask(ji,jj,1)                &
1361                        *(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1)    &
1362                        +frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))
1363                ENDDO
1364             ENDDO
1365                   
1366             ikchoix = 1
1367             CALL repcmo (frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)
1368             CALL lbc_lnk ('jpr_otx1', ztx2,'U', -1. )
1369             CALL lbc_lnk ('jpr_oty1', zty2,'V', -1. )
1370             frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)
1371             frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)
1372          ELSE
1373             CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
1374             frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1375             IF( srcv(jpr_otx2)%laction ) THEN
1376                CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
1377             ELSE
1378                CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
1379             ENDIF
1380          frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid 
1381               ENDIF
1382            ENDIF
1383            !                             
1384            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
1385               DO jj = 2, jpjm1                                          ! T ==> (U,V)
1386                  DO ji = fs_2, fs_jpim1   ! vector opt.
1387                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
1388                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
1389                  END DO
1390               END DO
1391               CALL lbc_lnk_multi( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1., frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
1392            ENDIF
1393            llnewtx = .TRUE.
1394         ELSE
1395            llnewtx = .FALSE.
1396         ENDIF
1397         !                                                   ! ========================= !
1398      ELSE                                                   !   No dynamical coupling   !
1399         !                                                   ! ========================= !
1400         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
1401         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
1402         llnewtx = .TRUE.
1403         !
1404      ENDIF
1405      !                                                      ! ========================= !
1406      !                                                      !    wind stress module     !   (taum)
1407      !                                                      ! ========================= !
1408      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
1409         ! => need to be done only when otx1 was changed
1410         IF( llnewtx ) THEN
1411            DO jj = 2, jpjm1
1412               DO ji = fs_2, fs_jpim1   ! vect. opt.
1413                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
1414                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
1415                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
1416               END DO
1417            END DO
1418            CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
1419            llnewtau = .TRUE.
1420         ELSE
1421            llnewtau = .FALSE.
1422         ENDIF
1423      ELSE
1424         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
1425         ! Stress module can be negative when received (interpolation problem)
1426         IF( llnewtau ) THEN
1427            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
1428         ENDIF
1429      ENDIF
1430      !
1431      !                                                      ! ========================= !
1432      !                                                      !      10 m wind speed      !   (wndm)
1433      !                                                      ! ========================= !
1434      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
1435         ! => need to be done only when taumod was changed
1436         IF( llnewtau ) THEN
1437            zcoef = 1. / ( zrhoa * zcdrag ) 
1438            DO jj = 1, jpj
1439               DO ji = 1, jpi 
1440                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
1441               END DO
1442            END DO
1443         ENDIF
1444      ENDIF
1445!!$      !                                                      ! ========================= !
1446!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      !
1447!!$      !                                                      ! ========================= !
1448!!$      cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1)
1449!!$      END SELECT
1450!!$
1451      zcloud_fra(:,:) = pp_cldf   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
1452      IF( ln_mixcpl ) THEN
1453         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:)
1454      ELSE
1455         cloud_fra(:,:) = zcloud_fra(:,:)
1456      ENDIF
1457      !                                                      ! ========================= !
1458      ! u(v)tau and taum will be modified by ice model
1459      ! -> need to be reset before each call of the ice/fsbc     
1460      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
1461         !
1462         IF( ln_mixcpl ) THEN
1463            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
1464            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
1465            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
1466            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
1467         ELSE
1468            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
1469            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
1470            taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
1471            wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
1472         ENDIF
1473         CALL iom_put( "taum_oce", taum )   ! output wind stress module
1474         
1475      ENDIF
1476
1477      IF (ln_medusa) THEN
1478        IF( srcv(jpr_atm_pco2)%laction) PCO2a_in_cpl(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
1479        IF( srcv(jpr_atm_dust)%laction) Dust_in_cpl(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
1480      ENDIF
1481      !
1482      !                                                      ! ================== !
1483      !                                                      ! atmosph. CO2 (ppm) !
1484      !                                                      ! ================== !
1485      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
1486      !
1487      !                                                      ! ========================= !
1488      !                                                      ! Mean Sea Level Pressure   !   (taum)
1489      !                                                      ! ========================= !
1490      IF( srcv(jpr_mslp)%laction ) THEN                    ! UKMO SHELF effect of atmospheric pressure on SSH
1491          IF( kt /= nit000 )   ssh_ibb(:,:) = ssh_ib(:,:)    !* Swap of ssh_ib fields
1492
1493          r1_grau = 1.e0 / (grav * rau0)               !* constant for optimization
1494          ssh_ib(:,:) = - ( frcv(jpr_mslp)%z3(:,:,1) - rpref ) * r1_grau    ! equivalent ssh (inverse barometer)
1495          apr   (:,:) =     frcv(jpr_mslp)%z3(:,:,1)                         !atmospheric pressure
1496   
1497          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)
1498      END IF 
1499      !
1500      IF( ln_sdw ) THEN  ! Stokes Drift correction activated
1501      !                                                      ! ========================= !
1502      !                                                      !       Stokes drift u      !
1503      !                                                      ! ========================= !
1504         IF( srcv(jpr_sdrftx)%laction ) ut0sd(:,:) = frcv(jpr_sdrftx)%z3(:,:,1)
1505      !
1506      !                                                      ! ========================= !
1507      !                                                      !       Stokes drift v      !
1508      !                                                      ! ========================= !
1509         IF( srcv(jpr_sdrfty)%laction ) vt0sd(:,:) = frcv(jpr_sdrfty)%z3(:,:,1)
1510      !
1511      !                                                      ! ========================= !
1512      !                                                      !      Wave mean period     !
1513      !                                                      ! ========================= !
1514         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)
1515      !
1516      !                                                      ! ========================= !
1517      !                                                      !  Significant wave height  !
1518      !                                                      ! ========================= !
1519         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
1520      !
1521      !                                                      ! ========================= ! 
1522      !                                                      !    Wave peak frequency    !
1523      !                                                      ! ========================= ! 
1524         IF( srcv(jpr_wfreq)%laction ) wfreq(:,:) = frcv(jpr_wfreq)%z3(:,:,1)
1525      !
1526      !                                                      ! ========================= !
1527      !                                                      !    Vertical mixing Qiao   !
1528      !                                                      ! ========================= !
1529         IF( srcv(jpr_wnum)%laction .AND. ln_zdfswm ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)
1530
1531         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
1532         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction &
1533                                      .OR. srcv(jpr_hsig)%laction   .OR. srcv(jpr_wfreq)%laction) THEN
1534            CALL sbc_stokes()
1535         ENDIF
1536      ENDIF
1537      !                                                      ! ========================= !
1538      !                                                      ! Stress adsorbed by waves  !
1539      !                                                      ! ========================= !
1540      IF( srcv(jpr_tauwoc)%laction .AND. ln_tauwoc ) tauoc_wave(:,:) = frcv(jpr_tauwoc)%z3(:,:,1)
1541
1542      !                                                      ! ========================= ! 
1543      !                                                      ! Stress component by waves !
1544      !                                                      ! ========================= ! 
1545      IF( srcv(jpr_tauwx)%laction .AND. srcv(jpr_tauwy)%laction .AND. ln_tauw ) THEN
1546         tauw_x(:,:) = frcv(jpr_tauwx)%z3(:,:,1)
1547         tauw_y(:,:) = frcv(jpr_tauwy)%z3(:,:,1)
1548      ENDIF
1549
1550      !                                                      ! ========================= !
1551      !                                                      !   Wave drag coefficient   !
1552      !                                                      ! ========================= !
1553      IF( srcv(jpr_wdrag)%laction .AND. ln_cdgw )   cdn_wave(:,:) = frcv(jpr_wdrag)%z3(:,:,1)
1554
1555      !  Fields received by SAS when OASIS coupling
1556      !  (arrays no more filled at sbcssm stage)
1557      !                                                      ! ================== !
1558      !                                                      !        SSS         !
1559      !                                                      ! ================== !
1560      IF( srcv(jpr_soce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1561         sss_m(:,:) = frcv(jpr_soce)%z3(:,:,1)
1562         CALL iom_put( 'sss_m', sss_m )
1563      ENDIF
1564      !                                               
1565      !                                                      ! ================== !
1566      !                                                      !        SST         !
1567      !                                                      ! ================== !
1568      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1569         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1)
1570         IF( srcv(jpr_soce)%laction .AND. l_useCT ) THEN    ! make sure that sst_m is the potential temperature
1571            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) )
1572         ENDIF
1573      ENDIF
1574      !                                                      ! ================== !
1575      !                                                      !        SSH         !
1576      !                                                      ! ================== !
1577      IF( srcv(jpr_ssh )%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1578         ssh_m(:,:) = frcv(jpr_ssh )%z3(:,:,1)
1579         CALL iom_put( 'ssh_m', ssh_m )
1580      ENDIF
1581      !                                                      ! ================== !
1582      !                                                      !  surface currents  !
1583      !                                                      ! ================== !
1584      IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1585         ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
1586         ub (:,:,1) = ssu_m(:,:)                             ! will be used in icestp in the call of ice_forcing_tau
1587         un (:,:,1) = ssu_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
1588         CALL iom_put( 'ssu_m', ssu_m )
1589      ENDIF
1590      IF( srcv(jpr_ocy1)%laction ) THEN
1591         ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
1592         vb (:,:,1) = ssv_m(:,:)                             ! will be used in icestp in the call of ice_forcing_tau
1593         vn (:,:,1) = ssv_m(:,:)                             ! will be used in sbc_cpl_snd if atmosphere coupling
1594         CALL iom_put( 'ssv_m', ssv_m )
1595      ENDIF
1596      !                                                      ! ======================== !
1597      !                                                      !  first T level thickness !
1598      !                                                      ! ======================== !
1599      IF( srcv(jpr_e3t1st )%laction ) THEN                   ! received by sas in case of opa <-> sas coupling
1600         e3t_m(:,:) = frcv(jpr_e3t1st )%z3(:,:,1)
1601         CALL iom_put( 'e3t_m', e3t_m(:,:) )
1602      ENDIF
1603      !                                                      ! ================================ !
1604      !                                                      !  fraction of solar net radiation !
1605      !                                                      ! ================================ !
1606      IF( srcv(jpr_fraqsr)%laction ) THEN                    ! received by sas in case of opa <-> sas coupling
1607         frq_m(:,:) = frcv(jpr_fraqsr)%z3(:,:,1)
1608         CALL iom_put( 'frq_m', frq_m )
1609      ENDIF
1610     
1611      !                                                      ! ========================= !
1612      IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
1613         !                                                   ! ========================= !
1614         !
1615         !                                                       ! total freshwater fluxes over the ocean (emp)
1616         IF( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction ) THEN
1617            SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
1618            CASE( 'conservative' )
1619               zemp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
1620            CASE( 'oce only', 'oce and ice' )
1621               zemp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
1622            CASE default
1623               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
1624            END SELECT
1625         ELSE
1626            zemp(:,:) = 0._wp
1627         ENDIF
1628         !
1629         !                                                        ! runoffs and calving (added in emp)
1630         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1631         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
1632 
1633         IF( srcv(jpr_icb)%laction )  THEN
1634             fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
1635             rnf(:,:)    = rnf(:,:) + fwficb(:,:)   ! iceberg added to runfofs
1636         ENDIF
1637         IF( srcv(jpr_isf)%laction )  fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)  ! fresh water flux from the isf (fwfisf <0 mean melting) 
1638       
1639         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
1640         ELSE                   ;   emp(:,:) =                              zemp(:,:)
1641         ENDIF
1642         !
1643         !                                                       ! non solar heat flux over the ocean (qns)
1644         IF(      srcv(jpr_qnsoce)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
1645         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
1646         ELSE                                       ;   zqns(:,:) = 0._wp
1647         END IF
1648         ! update qns over the free ocean with:
1649         IF( nn_components /= jp_iam_opa ) THEN
1650            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
1651            IF( srcv(jpr_snow  )%laction ) THEN
1652               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * rLfus   ! energy for melting solid precipitation over the free ocean
1653            ENDIF
1654         ENDIF
1655         !
1656         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus ! remove heat content associated to iceberg melting
1657         !
1658         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
1659         ELSE                   ;   qns(:,:) =                              zqns(:,:)
1660         ENDIF
1661
1662         !                                                       ! solar flux over the ocean          (qsr)
1663         IF     ( srcv(jpr_qsroce)%laction ) THEN   ;   zqsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
1664         ELSE IF( srcv(jpr_qsrmix)%laction ) then   ;   zqsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
1665         ELSE                                       ;   zqsr(:,:) = 0._wp
1666         ENDIF
1667         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
1668         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
1669         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
1670         ENDIF
1671         !
1672         ! salt flux over the ocean (received by opa in case of opa <-> sas coupling)
1673         IF( srcv(jpr_sflx )%laction )   sfx(:,:) = frcv(jpr_sflx  )%z3(:,:,1)
1674         ! Ice cover  (received by opa in case of opa <-> sas coupling)
1675         IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
1676         !
1677      ENDIF
1678
1679      !                                                        ! land ice masses : Greenland
1680      zepsilon = rn_iceshelf_fluxes_tolerance
1681
1682      IF( srcv(jpr_grnm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
1683     
1684         ! This is a zero dimensional, single value field.
1685         zgreenland_icesheet_mass_in =  frcv(jpr_grnm)%z3(1,1,1)
1686           
1687         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
1688
1689         IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
1690            ! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
1691            ! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
1692            zgreenland_icesheet_mass_b = zgreenland_icesheet_mass_in
1693            greenland_icesheet_mass = zgreenland_icesheet_mass_in
1694         ENDIF
1695
1696         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
1697            zgreenland_icesheet_mass_b = greenland_icesheet_mass
1698           
1699            ! Only update the mass if it has increased.
1700            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
1701               greenland_icesheet_mass = zgreenland_icesheet_mass_in
1702            ENDIF
1703           
1704            IF( zgreenland_icesheet_mass_b /= 0.0 ) &
1705           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed 
1706            greenland_icesheet_timelapsed = 0.0_wp       
1707         ENDIF
1708         IF(lwp .AND. ll_wrtstp) THEN
1709            WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
1710            WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass
1711            WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
1712            WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
1713         ENDIF
1714      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
1715         greenland_icesheet_mass_rate_of_change = rn_greenland_total_fw_flux
1716      ENDIF
1717
1718      !                                                        ! land ice masses : Antarctica
1719      IF( srcv(jpr_antm)%laction .AND. nn_coupled_iceshelf_fluxes == 1 ) THEN
1720         
1721         ! This is a zero dimensional, single value field.
1722         zantarctica_icesheet_mass_in = frcv(jpr_antm)%z3(1,1,1)
1723           
1724         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
1725
1726         IF( ln_iceshelf_init_atmos .AND. kt == 1 ) THEN
1727            ! On the first timestep (of an NRUN) force the ocean to ignore the icesheet masses in the ocean restart
1728            ! and take them from the atmosphere to avoid problems with using inconsistent ocean and atmosphere restarts.
1729            zantarctica_icesheet_mass_b = zantarctica_icesheet_mass_in
1730            antarctica_icesheet_mass = zantarctica_icesheet_mass_in
1731         ENDIF
1732
1733         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
1734            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
1735           
1736            ! Only update the mass if it has increased.
1737            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
1738               antarctica_icesheet_mass = zantarctica_icesheet_mass_in
1739            END IF
1740           
1741            IF( zantarctica_icesheet_mass_b /= 0.0 ) &
1742          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed 
1743            antarctica_icesheet_timelapsed = 0.0_wp       
1744         ENDIF
1745         IF(lwp .AND. ll_wrtstp) THEN
1746            WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
1747            WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass
1748            WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
1749            WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
1750         ENDIF
1751      ELSE IF ( nn_coupled_iceshelf_fluxes == 2 ) THEN
1752         antarctica_icesheet_mass_rate_of_change = rn_antarctica_total_fw_flux
1753      ENDIF
1754      !
1755      IF (ln_timing) CALL timing_stop('sbc_cpl_rcv')
1756   END SUBROUTINE sbc_cpl_rcv
1757   
1758
1759   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1760      !!----------------------------------------------------------------------
1761      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
1762      !!
1763      !! ** Purpose :   provide the stress over sea-ice in coupled mode
1764      !!
1765      !! ** Method  :   transform the received stress from the atmosphere into
1766      !!             an atmosphere-ice stress in the (i,j) ocean referencial
1767      !!             and at the velocity point of the sea-ice model:
1768      !!                'C'-grid : i- (j-) components given at U- (V-) point
1769      !!
1770      !!                The received stress are :
1771      !!                 - defined by 3 components (if cartesian coordinate)
1772      !!                        or by 2 components (if spherical)
1773      !!                 - oriented along geographical   coordinate (if eastward-northward)
1774      !!                        or  along the local grid coordinate (if local grid)
1775      !!                 - given at U- and V-point, resp.   if received on 2 grids
1776      !!                        or at a same point (T or I) if received on 1 grid
1777      !!                Therefore and if necessary, they are successively
1778      !!             processed in order to obtain them
1779      !!                 first  as  2 components on the sphere
1780      !!                 second as  2 components oriented along the local grid
1781      !!                 third  as  2 components on the ice grid point
1782      !!
1783      !!                Except in 'oce and ice' case, only one vector stress field
1784      !!             is received. It has already been processed in sbc_cpl_rcv
1785      !!             so that it is now defined as (i,j) components given at U-
1786      !!             and V-points, respectively. 
1787      !!
1788      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice
1789      !!----------------------------------------------------------------------
1790      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1791      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1792      !!
1793      INTEGER ::   ji, jj   ! dummy loop indices
1794      INTEGER ::   itx      ! index of taux over ice
1795      REAL(wp)                     ::   zztmp1, zztmp2
1796      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty 
1797      !!----------------------------------------------------------------------
1798      !
1799      IF (ln_timing) CALL timing_start('sbc_cpl_snd')
1800      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
1801      ELSE                                ;   itx =  jpr_otx1
1802      ENDIF
1803
1804      ! do something only if we just received the stress from atmosphere
1805      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
1806         !                                                      ! ======================= !
1807         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
1808            !                                                   ! ======================= !
1809           
1810            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1811               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1812               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
1813                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
1814               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1815               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1816               !
1817               IF( srcv(jpr_itx2)%laction ) THEN
1818                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
1819                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
1820                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1821                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1822               ENDIF
1823               !
1824            ENDIF
1825            !
1826            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1827               !                                                       ! (geographical to local grid -> rotate the components)
1828               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
1829               IF( srcv(jpr_itx2)%laction ) THEN
1830                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
1831               ELSE
1832                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
1833               ENDIF
1834               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1835               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
1836            ENDIF
1837            !                                                   ! ======================= !
1838         ELSE                                                   !     use ocean stress    !
1839            !                                                   ! ======================= !
1840            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
1841            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
1842            !
1843         ENDIF
1844         !                                                      ! ======================= !
1845         !                                                      !     put on ice grid     !
1846         !                                                      ! ======================= !
1847         !   
1848         !                                                  j+1   j     -----V---F
1849         ! ice stress on ice velocity point                              !       |
1850         ! (C-grid ==>(U,V))                                      j      |   T   U
1851         !                                                               |       |
1852         !                                                   j    j-1   -I-------|
1853         !                                               (for I)         |       |
1854         !                                                              i-1  i   i
1855         !                                                               i      i+1 (for I)
1856         SELECT CASE ( srcv(jpr_itx1)%clgrid )
1857         CASE( 'U' )
1858            p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1859            p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1860         CASE( 'T' )
1861            DO jj = 2, jpjm1                                   ! T ==> (U,V)
1862               DO ji = fs_2, fs_jpim1   ! vector opt.
1863                  ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology
1864                  zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
1865                  zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
1866                  p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1867                  p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1868               END DO
1869            END DO
1870            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
1871         END SELECT
1872         
1873      ENDIF
1874      !
1875   END SUBROUTINE sbc_cpl_ice_tau
1876   
1877
1878   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi )
1879      !!----------------------------------------------------------------------
1880      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1881      !!
1882      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
1883      !!
1884      !! ** Method  :   transform the fields received from the atmosphere into
1885      !!             surface heat and fresh water boundary condition for the
1886      !!             ice-ocean system. The following fields are provided:
1887      !!               * total non solar, solar and freshwater fluxes (qns_tot,
1888      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1889      !!             NB: emp_tot include runoffs and calving.
1890      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1891      !!             emp_ice = sublimation - solid precipitation as liquid
1892      !!             precipitation are re-routed directly to the ocean and
1893      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
1894      !!               * solid precipitation (sprecip), used to add to qns_tot
1895      !!             the heat lost associated to melting solid precipitation
1896      !!             over the ocean fraction.
1897      !!               * heat content of rain, snow and evap can also be provided,
1898      !!             otherwise heat flux associated with these mass flux are
1899      !!             guessed (qemp_oce, qemp_ice)
1900      !!
1901      !!             - the fluxes have been separated from the stress as
1902      !!               (a) they are updated at each ice time step compare to
1903      !!               an update at each coupled time step for the stress, and
1904      !!               (b) the conservative computation of the fluxes over the
1905      !!               sea-ice area requires the knowledge of the ice fraction
1906      !!               after the ice advection and before the ice thermodynamics,
1907      !!               so that the stress is updated before the ice dynamics
1908      !!               while the fluxes are updated after it.
1909      !!
1910      !! ** Details
1911      !!             qns_tot = (1-a) * qns_oce + a * qns_ice               => provided
1912      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
1913      !!
1914      !!             qsr_tot = (1-a) * qsr_oce + a * qsr_ice               => provided
1915      !!
1916      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce).
1917      !!                                                                      runoff (which includes rivers+icebergs) and iceshelf
1918      !!                                                                      are provided but not included in emp here. Only runoff will
1919      !!                                                                      be included in emp in other parts of NEMO code
1920      !! ** Action  :   update at each nf_ice time step:
1921      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1922      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1923      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
1924      !!                   emp_ice           ice sublimation - solid precipitation over the ice
1925      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
1926      !!                   sprecip           solid precipitation over the ocean 
1927      !!----------------------------------------------------------------------
1928      REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
1929      !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling
1930      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo
1931      REAL(wp), INTENT(in)   , DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
1932      REAL(wp), INTENT(inout), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] => inout for Met-Office
1933      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
1934      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
1935      !
1936      INTEGER  ::   ji, jj, jl   ! dummy loop index
1937      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
1938      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
1939      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
1940      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total
1941      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu
1942      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
1943      !!----------------------------------------------------------------------
1944      !
1945      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1946      ziceld(:,:) = 1._wp - picefr(:,:)
1947      zcptn (:,:) = rcp * sst_m(:,:)
1948
1949#if defined key_si3     
1950      !                                                      ! ========================= !
1951      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
1952      !                                                      ! ========================= !
1953      CASE ('coupled')
1954         IF (ln_scale_ice_flux) THEN
1955            WHERE( a_i(:,:,:) > 1.e-10_wp )
1956               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
1957               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
1958            ELSEWHERE
1959               qml_ice(:,:,:) = 0.0_wp
1960               qcn_ice(:,:,:) = 0.0_wp
1961            END WHERE
1962         ELSE
1963            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
1964            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
1965         ENDIF
1966      END SELECT
1967#endif
1968
1969      !
1970      !                                                      ! ========================= !
1971      !                                                      !    freshwater budget      !   (emp_tot)
1972      !                                                      ! ========================= !
1973      !
1974      !                                                           ! solid Precipitation                                (sprecip)
1975      !                                                           ! liquid + solid Precipitation                       (tprecip)
1976      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
1977      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
1978      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1979      CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1980         zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
1981         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
1982         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
1983      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1984         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1985         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
1986         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
1987         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
1988      CASE( 'none'      )       ! Not available as for now: needs additional coding below when computing zevap_oce
1989      !                         ! since fields received are not defined with none option
1990         CALL ctl_stop( 'STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
1991      END SELECT
1992
1993#if defined key_si3
1994
1995      ! --- evaporation over ice (kg/m2/s) --- !
1996      IF (ln_scale_ice_flux) THEN ! typically met-office requirements
1997         IF (sn_rcv_emp%clcat == 'yes') THEN
1998            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
1999            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp
2000            END WHERE
2001            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
2002            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
2003            END WHERE
2004         ELSE
2005            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:)
2006            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp
2007            END WHERE
2008            zevap_ice_total(:,:) = zevap_ice(:,:,1)
2009            DO jl = 2, jpl
2010               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
2011            ENDDO
2012         ENDIF
2013      ELSE
2014         IF (sn_rcv_emp%clcat == 'yes') THEN
2015            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
2016            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
2017            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
2018            END WHERE
2019         ELSE
2020            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1)
2021            zevap_ice_total(:,:) = zevap_ice(:,:,1)
2022            DO jl = 2, jpl
2023               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
2024            ENDDO
2025         ENDIF
2026      ENDIF
2027
2028      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
2029         ! For conservative case zemp_ice has not been defined yet. Do it now.
2030         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
2031      ENDIF
2032
2033      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
2034      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw )
2035     
2036      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
2037      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( picefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip
2038      zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice
2039
2040      ! --- evaporation over ocean (used later for qemp) --- !
2041      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
2042
2043      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
2044      ! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
2045      zdevap_ice(:,:) = 0._wp
2046     
2047      ! --- Continental fluxes --- !
2048      IF( srcv(jpr_rnf)%laction ) THEN   ! 2D runoffs (included in emp later on)
2049         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
2050      ENDIF
2051      IF( srcv(jpr_rnf_1d)%laction ) THEN ! 1D runoff
2052         CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:))
2053      ENDIF
2054      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot and emp_oce)
2055         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
2056         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
2057      ENDIF
2058      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
2059         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
2060         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
2061      ENDIF
2062      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
2063        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) 
2064      ENDIF
2065
2066      IF( ln_mixcpl ) THEN
2067         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
2068         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
2069         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:)
2070         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
2071         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
2072         DO jl = 1, jpl
2073            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:,jl) * zmsk(:,:)
2074            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:)    * zmsk(:,:)
2075         END DO
2076      ELSE
2077         emp_tot (:,:)   = zemp_tot (:,:)
2078         emp_ice (:,:)   = zemp_ice (:,:)
2079         emp_oce (:,:)   = zemp_oce (:,:)     
2080         sprecip (:,:)   = zsprecip (:,:)
2081         tprecip (:,:)   = ztprecip (:,:)
2082         evap_ice(:,:,:) = zevap_ice(:,:,:)
2083         DO jl = 1, jpl
2084            devap_ice(:,:,jl) = zdevap_ice(:,:)
2085         END DO
2086      ENDIF
2087
2088#else
2089      zsnw(:,:) = picefr(:,:)
2090      ! --- Continental fluxes --- !
2091      IF( srcv(jpr_rnf)%laction ) THEN   ! 2D runoffs (included in emp later on)
2092         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
2093      ENDIF
2094      IF( srcv(jpr_rnf_1d)%laction ) THEN  ! 1D runoff
2095         CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:)) 
2096      ENDIF
2097      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot)
2098         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
2099      ENDIF
2100      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
2101         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
2102         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
2103      ENDIF
2104      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf <0 mean melting)
2105        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1)
2106      ENDIF
2107      !
2108      IF( ln_mixcpl ) THEN
2109         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
2110         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
2111         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
2112         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
2113      ELSE
2114         emp_tot(:,:) =                                  zemp_tot(:,:)
2115         emp_ice(:,:) =                                  zemp_ice(:,:)
2116         sprecip(:,:) =                                  zsprecip(:,:)
2117         tprecip(:,:) =                                  ztprecip(:,:)
2118      ENDIF
2119      !
2120#endif
2121
2122      ! outputs
2123!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
2124!!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
2125      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
2126      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
2127      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
2128      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
2129      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation
2130      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
2131      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
2132      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average)
2133      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )     ! Sublimation over sea-ice (cell average)
2134      IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
2135         &                                                         - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average)
2136      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
2137      !
2138      !                                                      ! ========================= !
2139      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
2140      !                                                      ! ========================= !
2141      CASE( 'oce only' )         ! the required field is directly provided
2142
2143         ! Get the sea ice non solar heat flux from conductive and melting fluxes
2144         IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
2145            zqns_ice(:,:,:) = qml_ice(:,:,:) + qcn_ice(:,:,:)
2146         ELSE
2147            zqns_ice(:,:,:) = 0._wp
2148         ENDIF
2149
2150         ! Calculate the total non solar heat flux. The ocean only non solar heat flux (zqns_oce) will be recalculated after this CASE
2151         ! statement to be consistent with other coupling methods even though .zqns_oce = frcv(jpr_qnsoce)%z3(:,:,1)
2152         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) + SUM( zqns_ice(:,:,:) * a_i(:,:,:), dim=3 )
2153         
2154      CASE( 'conservative' )     ! the required fields are directly provided
2155         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
2156         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
2157            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
2158         ELSE
2159            DO jl = 1, jpl
2160               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
2161            END DO
2162         ENDIF
2163      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
2164         zqns_tot(:,:) =  ziceld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
2165         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
2166            DO jl=1,jpl
2167               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
2168               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
2169            ENDDO
2170         ELSE
2171            zqns_tot(:,:) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
2172            DO jl = 1, jpl
2173               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
2174            END DO
2175         ENDIF
2176      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
2177! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
2178         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
2179         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
2180            DO jl = 1, jpl
2181               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    &
2182                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
2183                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
2184            END DO
2185         ELSE
2186            DO jl = 1, jpl
2187               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    &
2188                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
2189                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
2190            END DO
2191         ENDIF
2192      END SELECT
2193      !                                     
2194      ! --- calving (removed from qns_tot) --- !
2195      IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * rLfus  ! remove latent heat of calving
2196                                                                                                     ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
2197      ! --- iceberg (removed from qns_tot) --- !
2198      IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus  ! remove latent heat of iceberg melting
2199
2200#if defined key_si3     
2201      ! --- non solar flux over ocean --- !
2202      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
2203      zqns_oce = 0._wp
2204      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
2205
2206      ! Heat content per unit mass of snow (J/kg)
2207      WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
2208      ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:)
2209      ENDWHERE
2210      ! Heat content per unit mass of rain (J/kg)
2211      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * ziceld(:,:) ) 
2212
2213      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
2214      zqprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus )
2215
2216      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
2217      DO jl = 1, jpl
2218         ! zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account
2219
2220         ! How much enthalpy is stored in sublimating snow and ice. At this stage we don't know if it is snow or ice that is
2221         ! sublimating so we will use the combined snow and ice layer temperature t1_ice.
2222         zqevap_ice(:,:,jl) = -zevap_ice(:,:,jl) * ( ( rt0 - t1_ice(:,:,jl) ) * rcpi + rLfus )
2223
2224      END DO
2225
2226      ! --- heat flux associated with emp (W/m2) --- !
2227      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn   (:,:)   &        ! evap
2228         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptrain(:,:)   &        ! liquid precip
2229         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus )  ! solid precip over ocean + snow melting
2230      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - rLfus )  ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
2231!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
2232!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhos  ! solid precip over ice
2233     
2234      ! --- total non solar flux (including evap/precip) --- !
2235      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)
2236
2237      ! --- in case both coupled/forced are active, we must mix values --- !
2238      IF( ln_mixcpl ) THEN
2239         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
2240         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
2241         DO jl=1,jpl
2242            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
2243            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
2244         ENDDO
2245         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
2246         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
2247         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
2248      ELSE
2249         qns_tot  (:,:  ) = zqns_tot  (:,:  )
2250         qns_oce  (:,:  ) = zqns_oce  (:,:  )
2251         qns_ice  (:,:,:) = zqns_ice  (:,:,:)
2252         qevap_ice(:,:,:) = zqevap_ice(:,:,:)
2253         qprec_ice(:,:  ) = zqprec_ice(:,:  )
2254         qemp_oce (:,:  ) = zqemp_oce (:,:  )
2255         qemp_ice (:,:  ) = zqemp_ice (:,:  )
2256      ENDIF
2257
2258#else
2259      zcptsnw (:,:) = zcptn(:,:)
2260      zcptrain(:,:) = zcptn(:,:)
2261     
2262      ! clem: this formulation is certainly wrong... but better than it was...
2263      zqns_tot(:,:) = zqns_tot(:,:)                             &          ! zqns_tot update over free ocean with:
2264         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &          ! remove the latent heat flux of solid precip. melting
2265         &          - (  zemp_tot(:,:)                          &          ! remove the heat content of mass flux (assumed to be at SST)
2266         &             - zemp_ice(:,:) ) * zcptn(:,:) 
2267
2268     IF( ln_mixcpl ) THEN
2269         qns_tot(:,:) = qns(:,:) * ziceld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
2270         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
2271         DO jl=1,jpl
2272            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
2273         ENDDO
2274      ELSE
2275         qns_tot(:,:  ) = zqns_tot(:,:  )
2276         qns_ice(:,:,:) = zqns_ice(:,:,:)
2277      ENDIF
2278
2279#endif
2280      ! outputs
2281      IF ( srcv(jpr_cal)%laction ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )   ! latent heat from calving
2282      IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )   ! latent heat from icebergs melting
2283      IF ( iom_use(   'hflx_rain_cea') ) CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average)
2284      IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) )  &
2285           &                         * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
2286      IF ( iom_use(   'hflx_prec_cea') ) CALL iom_put('hflx_prec_cea' ,    sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  &                    ! heat flux from all precip (cell avg)
2287         &                          + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
2288      IF ( iom_use(   'hflx_snow_cea') ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average)
2289      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) )   ! heat flux from snow (over ocean)
2290      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *  zsnw(:,:) )              ! heat flux from snow (over ice)
2291      ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
2292      !
2293      !                                                      ! ========================= !
2294      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
2295      !                                                      ! ========================= !
2296      CASE( 'oce only' )
2297         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
2298
2299         ! For the Met Office the only sea ice solar flux is the transmitted qsr which is added onto zqsr_ice
2300         ! further down. Therefore start zqsr_ice off at zero.
2301         zqsr_ice(:,:,:) = 0._wp
2302      CASE( 'conservative' )
2303         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
2304         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
2305            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
2306         ELSE
2307            ! Set all category values equal for the moment
2308            DO jl = 1, jpl
2309               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
2310            END DO
2311         ENDIF
2312      CASE( 'oce and ice' )
2313         zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
2314         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
2315            DO jl = 1, jpl
2316               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
2317               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
2318            END DO
2319         ELSE
2320            zqsr_tot(:,:) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
2321            DO jl = 1, jpl
2322               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
2323            END DO
2324         ENDIF
2325      CASE( 'mixed oce-ice' )
2326         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
2327! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
2328!       Create solar heat flux over ice using incoming solar heat flux and albedos
2329!       ( see OASIS3 user guide, 5th edition, p39 )
2330         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
2331            DO jl = 1, jpl
2332               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:,jl) * ( 1.- palbi(:,:,jl) )   &
2333                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
2334                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
2335            END DO
2336         ELSE
2337            DO jl = 1, jpl
2338               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:, 1) * ( 1.- palbi(:,:,jl) )   &
2339                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
2340                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
2341            END DO
2342         ENDIF
2343      CASE( 'none'      )       ! Not available as for now: needs additional coding 
2344      !                         ! since fields received, here zqsr_tot,  are not defined with none option
2345         CALL ctl_stop( 'STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_qsr value in namelist namsbc_cpl' )
2346      END SELECT
2347      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
2348         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
2349         DO jl = 1, jpl
2350            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
2351         END DO
2352      ENDIF
2353
2354#if defined key_si3     
2355      !                                                      ! ========================= !
2356      !                                                      !      Transmitted Qsr      !   [W/m2]
2357      !                                                      ! ========================= !
2358      IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==!
2359         !
2360         IF( nn_qtrice == 0 ) THEN
2361            ! formulation derived from Grenfell and Maykut (1977), where transmission rate
2362            !    1) depends on cloudiness
2363            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
2364            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
2365            !    2) is 0 when there is any snow
2366            !    3) tends to 1 for thin ice
2367            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
2368            DO jl = 1, jpl
2369               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
2370                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
2371               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm
2372                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:)
2373               ELSEWHERE                                                           ! zero when hs>0
2374                  zqtr_ice_top(:,:,jl) = 0._wp 
2375               END WHERE
2376            ENDDO
2377         ELSEIF( nn_qtrice == 1 ) THEN
2378            ! formulation is derived from the thesis of M. Lebrun (2019).
2379            !    It represents the best fit using several sets of observations
2380            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
2381            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:)
2382         ENDIF
2383         !     
2384      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
2385         !
2386         SELECT CASE( TRIM( sn_rcv_qtr%cldes ) )
2387         !
2388         !      ! ===> here we receive the qtr_ice_top array from the coupler
2389         CASE ('coupled')
2390            IF (ln_scale_ice_flux) THEN
2391               WHERE( a_i(:,:,:) > 0.0_wp ) zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:)) * a_i_last_couple(:,:,:) / a_i(:,:,:)
2392               WHERE( a_i(:,:,:) <= 0.0_wp ) zqtr_ice_top(:,:,:) = 0.0_wp
2393            ELSE
2394               zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:))
2395            ENDIF
2396
2397            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
2398            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
2399            zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )
2400
2401         !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
2402         CASE ('none')
2403            zqtr_ice_top(:,:,:) = 0._wp
2404         END SELECT
2405
2406         !
2407      ENDIF
2408      !
2409      IF( ln_mixcpl ) THEN
2410         DO jl=1,jpl
2411            qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:)
2412         ENDDO
2413      ELSE
2414         qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:)
2415      ENDIF
2416
2417      ! --- solar flux over ocean --- !
2418      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
2419      zqsr_oce = 0._wp
2420      WHERE( ziceld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / ziceld(:,:)
2421
2422      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
2423      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
2424#endif
2425
2426      IF( ln_mixcpl ) THEN
2427         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
2428         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
2429         DO jl = 1, jpl
2430            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
2431         END DO
2432      ELSE
2433         qsr_tot(:,:  ) = zqsr_tot(:,:  )
2434         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
2435      ENDIF
2436
2437      !                                                      ! ========================= !
2438      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
2439      !                                                      ! ========================= !
2440      CASE ('coupled')
2441         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
2442            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
2443         ELSE
2444            ! Set all category values equal for the moment
2445            DO jl=1,jpl
2446               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
2447            ENDDO
2448         ENDIF
2449      CASE( 'none' ) 
2450         zdqns_ice(:,:,:) = 0._wp
2451      END SELECT
2452     
2453      IF( ln_mixcpl ) THEN
2454         DO jl=1,jpl
2455            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
2456         ENDDO
2457      ELSE
2458         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
2459      ENDIF
2460
2461#if defined key_si3
2462      !                                                      ! ================== !
2463      !                                                      !   ice skin temp.   !
2464      !                                                      ! ================== !
2465      ! needed by Met Office
2466      IF( srcv(jpr_ts_ice)%laction ) THEN
2467         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   ztsu(:,:,:) =   0. + rt0 
2468         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   ztsu(:,:,:) = -60. + rt0
2469         ELSEWHERE                                        ;   ztsu(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) + rt0
2470         END WHERE
2471         !
2472         IF( ln_mixcpl ) THEN
2473            DO jl=1,jpl
2474               pist(:,:,jl) = pist(:,:,jl) * xcplmask(:,:,0) + ztsu(:,:,jl) * zmsk(:,:)
2475            ENDDO
2476         ELSE
2477            pist(:,:,:) = ztsu(:,:,:)
2478         ENDIF
2479         !
2480      ENDIF
2481      !
2482#endif
2483      !
2484   END SUBROUTINE sbc_cpl_ice_flx
2485   
2486   
2487   SUBROUTINE sbc_cpl_snd( kt )
2488      !!----------------------------------------------------------------------
2489      !!             ***  ROUTINE sbc_cpl_snd  ***
2490      !!
2491      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
2492      !!
2493      !! ** Method  :   send to the atmosphere through a call to cpl_snd
2494      !!              all the needed fields (as defined in sbc_cpl_init)
2495      !!----------------------------------------------------------------------
2496      INTEGER, INTENT(in) ::   kt
2497      !
2498      INTEGER ::   ji, jj, jl   ! dummy loop indices
2499      INTEGER ::   ikchoix
2500      INTEGER ::   isec, info   ! local integer
2501      REAL(wp) ::   zumax, zvmax
2502      REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
2503      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4   
2504      !!----------------------------------------------------------------------
2505      !
2506      isec = ( kt - nit000 ) * NINT( rdt )        ! date of exchanges
2507
2508      zfr_l(:,:) = 1.- fr_i(:,:)
2509      !                                                      ! ------------------------- !
2510      !                                                      !    Surface temperature    !   in Kelvin
2511      !                                                      ! ------------------------- !
2512      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
2513         
2514         IF ( nn_components == jp_iam_opa ) THEN
2515            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
2516         ELSE
2517            ! we must send the surface potential temperature
2518            IF( l_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
2519            ELSE                   ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
2520            ENDIF
2521            !
2522            SELECT CASE( sn_snd_temp%cldes)
2523            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2524            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2525               SELECT CASE( sn_snd_temp%clcat )
2526               CASE( 'yes' )   
2527                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
2528               CASE( 'no' )
2529                  WHERE( SUM( a_i, dim=3 ) /= 0. )
2530                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
2531                  ELSEWHERE
2532                     ztmp3(:,:,1) = rt0
2533                  END WHERE
2534               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2535               END SELECT
2536            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
2537               SELECT CASE( sn_snd_temp%clcat )
2538               CASE( 'yes' )   
2539                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2540               CASE( 'no' )
2541                  ztmp3(:,:,:) = 0.0
2542                  DO jl=1,jpl
2543                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
2544                  ENDDO
2545               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2546               END SELECT
2547            CASE( 'oce and weighted ice')    ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0 
2548               SELECT CASE( sn_snd_temp%clcat ) 
2549               CASE( 'yes' )   
2550                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
2551               CASE( 'no' ) 
2552                  ztmp3(:,:,:) = 0.0 
2553                  DO jl=1,jpl 
2554                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl) 
2555                  ENDDO 
2556               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 
2557               END SELECT
2558            CASE( 'mixed oce-ice'        )   
2559               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
2560               DO jl=1,jpl
2561                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
2562               ENDDO
2563            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
2564            END SELECT
2565         ENDIF
2566         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2567         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
2568         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2569      ENDIF
2570      !
2571      !                                                      ! ------------------------- !
2572      !                                                      ! 1st layer ice/snow temp.  !
2573      !                                                      ! ------------------------- !
2574#if defined key_si3
2575      ! needed by  Met Office
2576      IF( ssnd(jps_ttilyr)%laction) THEN
2577         SELECT CASE( sn_snd_ttilyr%cldes)
2578         CASE ('weighted ice')
2579            ztmp3(:,:,1:jpl) = t1_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
2580         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ttilyr%cldes' )
2581         END SELECT
2582         IF( ssnd(jps_ttilyr)%laction )   CALL cpl_snd( jps_ttilyr, isec, ztmp3, info )
2583      ENDIF
2584#endif
2585      !                                                      ! ------------------------- !
2586      !                                                      !           Albedo          !
2587      !                                                      ! ------------------------- !
2588      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
2589          SELECT CASE( sn_snd_alb%cldes )
2590          CASE( 'ice' )
2591             SELECT CASE( sn_snd_alb%clcat )
2592             CASE( 'yes' )   
2593                ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
2594             CASE( 'no' )
2595                WHERE( SUM( a_i, dim=3 ) /= 0. )
2596                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
2597                ELSEWHERE
2598                   ztmp1(:,:) = alb_oce_mix(:,:)
2599                END WHERE
2600             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
2601             END SELECT
2602          CASE( 'weighted ice' )   ;
2603             SELECT CASE( sn_snd_alb%clcat )
2604             CASE( 'yes' )   
2605                ztmp3(:,:,1:jpl) =  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2606             CASE( 'no' )
2607                WHERE( fr_i (:,:) > 0. )
2608                   ztmp1(:,:) = SUM (  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
2609                ELSEWHERE
2610                   ztmp1(:,:) = 0.
2611                END WHERE
2612             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ice%clcat' )
2613             END SELECT
2614          CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
2615         END SELECT
2616
2617         SELECT CASE( sn_snd_alb%clcat )
2618            CASE( 'yes' )   
2619               CALL cpl_snd( jps_albice, isec, ztmp3, info )      !-> MV this has never been checked in coupled mode
2620            CASE( 'no'  )   
2621               CALL cpl_snd( jps_albice, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
2622         END SELECT
2623      ENDIF
2624
2625      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
2626         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
2627         DO jl = 1, jpl
2628            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
2629         END DO
2630         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2631      ENDIF
2632      !                                                      ! ------------------------- !
2633      !                                                      !  Ice fraction & Thickness !
2634      !                                                      ! ------------------------- !
2635      ! Send ice fraction field to atmosphere
2636      IF( ssnd(jps_fice)%laction ) THEN
2637         SELECT CASE( sn_snd_thick%clcat )
2638         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
2639         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
2640         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2641         END SELECT
2642         CALL cpl_snd( jps_fice, isec, ztmp3, info )
2643      ENDIF
2644
2645#if defined key_si3 || defined key_cice
2646      ! If this coupling was successful then save ice fraction for use between coupling points.
2647      ! This is needed for some calculations where the ice fraction at the last coupling point
2648      ! is needed.
2649      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. & 
2650         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN
2651         IF ( sn_snd_thick%clcat == 'yes' ) THEN
2652           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl)
2653         ENDIF
2654      ENDIF
2655#endif
2656
2657      IF( ssnd(jps_fice1)%laction ) THEN
2658         SELECT CASE( sn_snd_thick1%clcat )
2659         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
2660         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
2661         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' )
2662         END SELECT
2663         CALL cpl_snd( jps_fice1, isec, ztmp3, info )
2664      ENDIF
2665     
2666      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
2667      IF( ssnd(jps_fice2)%laction ) THEN
2668         ztmp3(:,:,1) = fr_i(:,:)
2669         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
2670      ENDIF
2671
2672      ! Send ice and snow thickness field
2673      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
2674         SELECT CASE( sn_snd_thick%cldes)
2675         CASE( 'none'                  )       ! nothing to do
2676         CASE( 'weighted ice and snow' )   
2677            SELECT CASE( sn_snd_thick%clcat )
2678            CASE( 'yes' )   
2679               ztmp3(:,:,1:jpl) =  h_i(:,:,1:jpl) * a_i(:,:,1:jpl)
2680               ztmp4(:,:,1:jpl) =  h_s(:,:,1:jpl) * a_i(:,:,1:jpl)
2681            CASE( 'no' )
2682               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
2683               DO jl=1,jpl
2684                  ztmp3(:,:,1) = ztmp3(:,:,1) + h_i(:,:,jl) * a_i(:,:,jl)
2685                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_s(:,:,jl) * a_i(:,:,jl)
2686               ENDDO
2687            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2688            END SELECT
2689         CASE( 'ice and snow'         )   
2690            SELECT CASE( sn_snd_thick%clcat )
2691            CASE( 'yes' )
2692               ztmp3(:,:,1:jpl) = h_i(:,:,1:jpl)
2693               ztmp4(:,:,1:jpl) = h_s(:,:,1:jpl)
2694            CASE( 'no' )
2695               WHERE( SUM( a_i, dim=3 ) /= 0. )
2696                  ztmp3(:,:,1) = SUM( h_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
2697                  ztmp4(:,:,1) = SUM( h_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
2698               ELSEWHERE
2699                 ztmp3(:,:,1) = 0.
2700                 ztmp4(:,:,1) = 0.
2701               END WHERE
2702            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2703            END SELECT
2704         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
2705         END SELECT
2706         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
2707         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
2708      ENDIF
2709
2710#if defined key_si3
2711      !                                                      ! ------------------------- !
2712      !                                                      !      Ice melt ponds       !
2713      !                                                      ! ------------------------- !
2714      ! needed by Met Office: 1) fraction of ponded ice 2) local/actual pond depth
2715      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN
2716         SELECT CASE( sn_snd_mpnd%cldes) 
2717         CASE( 'ice only' ) 
2718            SELECT CASE( sn_snd_mpnd%clcat ) 
2719            CASE( 'yes' ) 
2720               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl)
2721               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl) 
2722            CASE( 'no' ) 
2723               ztmp3(:,:,:) = 0.0 
2724               ztmp4(:,:,:) = 0.0 
2725               DO jl=1,jpl 
2726                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)
2727                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)
2728               ENDDO 
2729            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' ) 
2730            END SELECT 
2731         CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%cldes' )     
2732         END SELECT 
2733         IF( ssnd(jps_a_p)%laction  )   CALL cpl_snd( jps_a_p , isec, ztmp3, info )     
2734         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )     
2735      ENDIF 
2736      !
2737      !                                                      ! ------------------------- !
2738      !                                                      !     Ice conductivity      !
2739      !                                                      ! ------------------------- !
2740      ! needed by Met Office
2741      IF( ssnd(jps_kice)%laction ) THEN
2742         SELECT CASE( sn_snd_cond%cldes) 
2743         CASE( 'weighted ice' )   
2744            SELECT CASE( sn_snd_cond%clcat ) 
2745            CASE( 'yes' )   
2746          ztmp3(:,:,1:jpl) =  cnd_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
2747            CASE( 'no' ) 
2748               ztmp3(:,:,:) = 0.0 
2749               DO jl=1,jpl 
2750                 ztmp3(:,:,1) = ztmp3(:,:,1) + cnd_ice(:,:,jl) * a_i(:,:,jl) 
2751               ENDDO 
2752            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' ) 
2753            END SELECT
2754         CASE( 'ice only' )   
2755           ztmp3(:,:,1:jpl) = cnd_ice(:,:,1:jpl) 
2756         CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%cldes' )     
2757         END SELECT
2758         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info ) 
2759      ENDIF 
2760#endif
2761
2762      !                                                      ! ------------------------- !
2763      !                                                      !  CO2 flux from PISCES     !
2764      !                                                      ! ------------------------- !
2765      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN 
2766         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s
2767         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info ) 
2768      ENDIF 
2769      !
2770      IF (ln_medusa) THEN
2771      !                                                      ! ---------------------------------------------- !
2772      !                                                      !  CO2 flux, DMS and chlorophyll from MEDUSA     !
2773      !                                                      ! ---------------------------------------------- !
2774         IF ( ssnd(jps_bio_co2)%laction ) THEN
2775            CALL cpl_snd( jps_bio_co2, isec, RESHAPE( CO2Flux_out_cpl, (/jpi,jpj,1/) ), info )
2776         ENDIF
2777
2778         IF ( ssnd(jps_bio_dms)%laction )  THEN
2779            CALL cpl_snd( jps_bio_dms, isec, RESHAPE( DMS_out_cpl, (/jpi,jpj,1/) ), info )
2780         ENDIF
2781
2782         IF ( ssnd(jps_bio_chloro)%laction )  THEN
2783            CALL cpl_snd( jps_bio_chloro, isec, RESHAPE( chloro_out_cpl, (/jpi,jpj,1/) ), info )
2784         ENDIF
2785      ENDIF
2786      !
2787      !                                                      ! ------------------------- !
2788      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
2789         !                                                   ! ------------------------- !
2790         !   
2791         !                                                  j+1   j     -----V---F
2792         ! surface velocity always sent from T point                     !       |
2793         ! [except for HadGEM3]                                   j      |   T   U
2794         !                                                               |       |
2795         !                                                   j    j-1   -I-------|
2796         !                                               (for I)         |       |
2797         !                                                              i-1  i   i
2798         !                                                               i      i+1 (for I)
2799         IF( nn_components == jp_iam_opa ) THEN
2800            zotx1(:,:) = un(:,:,1) 
2801            zoty1(:,:) = vn(:,:,1) 
2802         ELSE       
2803            SELECT CASE( TRIM( sn_snd_crt%cldes ) )
2804            CASE( 'oce only'             )      ! C-grid ==> T
2805               IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
2806                  DO jj = 2, jpjm1
2807                     DO ji = fs_2, fs_jpim1   ! vector opt.
2808                        zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
2809                        zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
2810                     END DO
2811                  END DO
2812               ELSE
2813! Temporarily Changed for UKV
2814                  DO jj = 2, jpjm1
2815                     DO ji = 2, jpim1
2816                        zotx1(ji,jj) = un(ji,jj,1)
2817                        zoty1(ji,jj) = vn(ji,jj,1)
2818                     END DO
2819                  END DO
2820               ENDIF
2821            CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T 
2822               DO jj = 2, jpjm1
2823                  DO ji = fs_2, fs_jpim1   ! vector opt.
2824                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2825                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
2826                     zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2827                     zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2828                  END DO
2829               END DO
2830               CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1., zity1, 'T', -1. )
2831            CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T
2832               DO jj = 2, jpjm1
2833                  DO ji = fs_2, fs_jpim1   ! vector opt.
2834                     zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
2835                        &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2836                     zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
2837                        &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2838                  END DO
2839               END DO
2840            END SELECT
2841            CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.,  zoty1, ssnd(jps_ocy1)%clgrid, -1. )
2842            !
2843         ENDIF
2844         !
2845         !
2846         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2847            !                                                                     ! Ocean component
2848            IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
2849               CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
2850               CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
2851               zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
2852               zoty1(:,:) = ztmp2(:,:)
2853               IF( ssnd(jps_ivx1)%laction ) THEN                                  ! Ice component
2854                  CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
2855                  CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
2856                  zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
2857                  zity1(:,:) = ztmp2(:,:)
2858               ENDIF
2859            ELSE
2860               ! Temporary code for HadGEM3 - will be removed eventually.
2861               ! Only applies when we want uvel on U grid and vvel on V grid
2862               ! Rotate U and V onto geographic grid before sending.
2863
2864               DO jj=2,jpjm1
2865                  DO ji=2,jpim1
2866                     ztmp1(ji,jj)=0.25*vmask(ji,jj,1)                  &
2867                          *(zotx1(ji,jj)+zotx1(ji-1,jj)    &
2868                          +zotx1(ji,jj+1)+zotx1(ji-1,jj+1))
2869                     ztmp2(ji,jj)=0.25*umask(ji,jj,1)                  &
2870                          *(zoty1(ji,jj)+zoty1(ji+1,jj)    &
2871                          +zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
2872                  ENDDO
2873               ENDDO
2874               
2875               ! Ensure any N fold and wrap columns are updated
2876               CALL lbc_lnk('zotx1', ztmp1, 'V', -1.0)
2877               CALL lbc_lnk('zoty1', ztmp2, 'U', -1.0)
2878               
2879               ikchoix = -1
2880               CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix)
2881           ENDIF
2882         ENDIF
2883         !
2884         ! spherical coordinates to cartesian -> 2 components to 3 components
2885         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
2886            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2887            ztmp2(:,:) = zoty1(:,:)
2888            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2889            !
2890            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2891               ztmp1(:,:) = zitx1(:,:)
2892               ztmp1(:,:) = zity1(:,:)
2893               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2894            ENDIF
2895         ENDIF
2896         !
2897         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2898         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2899         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
2900         !
2901         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
2902         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
2903         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
2904         !
2905      ENDIF
2906      !
2907      !                                                      ! ------------------------- !
2908      !                                                      !  Surface current to waves !
2909      !                                                      ! ------------------------- !
2910      IF( ssnd(jps_ocxw)%laction .OR. ssnd(jps_ocyw)%laction ) THEN 
2911          !     
2912          !                                                  j+1  j     -----V---F
2913          ! surface velocity always sent from T point                    !       |
2914          !                                                       j      |   T   U
2915          !                                                              |       |
2916          !                                                   j   j-1   -I-------|
2917          !                                               (for I)        |       |
2918          !                                                             i-1  i   i
2919          !                                                              i      i+1 (for I)
2920          SELECT CASE( TRIM( sn_snd_crtw%cldes ) ) 
2921          CASE( 'oce only'             )      ! C-grid ==> T
2922             DO jj = 2, jpjm1 
2923                DO ji = fs_2, fs_jpim1   ! vector opt.
2924                   zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) ) 
2925                   zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji , jj-1,1) ) 
2926                END DO
2927             END DO
2928          CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T   
2929             DO jj = 2, jpjm1 
2930                DO ji = fs_2, fs_jpim1   ! vector opt.
2931                   zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   
2932                   zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2933                   zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
2934                   zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
2935                END DO
2936             END DO
2937             CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.,  zity1, 'T', -1. ) 
2938          CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T 
2939             DO jj = 2, jpjm1 
2940                DO ji = fs_2, fs_jpim1   ! vector opt.
2941                   zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   & 
2942                      &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj) 
2943                   zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2944                      &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj) 
2945                END DO
2946             END DO
2947          END SELECT
2948         CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1., zoty1, ssnd(jps_ocyw)%clgrid, -1. ) 
2949         !
2950         !
2951         IF( TRIM( sn_snd_crtw%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2952         !                                                                        ! Ocean component
2953            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->e', ztmp1 )       ! 1st component 
2954            CALL rot_rep( zotx1, zoty1, ssnd(jps_ocxw)%clgrid, 'ij->n', ztmp2 )       ! 2nd component 
2955            zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components 
2956            zoty1(:,:) = ztmp2(:,:) 
2957            IF( ssnd(jps_ivx1)%laction ) THEN                                     ! Ice component
2958               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component 
2959               CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component 
2960               zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components 
2961               zity1(:,:) = ztmp2(:,:) 
2962            ENDIF
2963         ENDIF 
2964         !
2965!         ! spherical coordinates to cartesian -> 2 components to 3 components
2966!         IF( TRIM( sn_snd_crtw%clvref ) == 'cartesian' ) THEN
2967!            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2968!            ztmp2(:,:) = zoty1(:,:)
2969!            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2970!            !
2971!            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2972!               ztmp1(:,:) = zitx1(:,:)
2973!               ztmp1(:,:) = zity1(:,:)
2974!               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2975!            ENDIF
2976!         ENDIF
2977         !
2978         IF( ssnd(jps_ocxw)%laction )   CALL cpl_snd( jps_ocxw, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2979         IF( ssnd(jps_ocyw)%laction )   CALL cpl_snd( jps_ocyw, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2980         
2981      ENDIF 
2982      !
2983      IF( ssnd(jps_ficet)%laction ) THEN
2984         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info ) 
2985      END IF 
2986      !                                                      ! ------------------------- !
2987      !                                                      !   Water levels to waves   !
2988      !                                                      ! ------------------------- !
2989      IF( ssnd(jps_wlev)%laction ) THEN
2990         IF( ln_apr_dyn ) THEN 
2991            IF( kt /= nit000 ) THEN 
2992               ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
2993            ELSE 
2994               ztmp1(:,:) = sshb(:,:) 
2995            ENDIF 
2996         ELSE 
2997            ztmp1(:,:) = sshn(:,:) 
2998         ENDIF 
2999         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
3000      END IF 
3001      !
3002      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
3003      !                                                        ! SSH
3004      IF( ssnd(jps_ssh )%laction )  THEN
3005         !                          ! removed inverse barometer ssh when Patm
3006         !                          forcing is used (for sea-ice dynamics)
3007         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
3008         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
3009         ENDIF
3010         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
3011
3012      ENDIF
3013      !                                                        ! SSS
3014      IF( ssnd(jps_soce  )%laction )  THEN
3015         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
3016      ENDIF
3017      !                                                        ! first T level thickness
3018      IF( ssnd(jps_e3t1st )%laction )  THEN
3019         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( e3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
3020      ENDIF
3021      !                                                        ! Qsr fraction
3022      IF( ssnd(jps_fraqsr)%laction )  THEN
3023         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
3024      ENDIF
3025      !
3026      !  Fields sent by SAS to OPA when OASIS coupling
3027      !                                                        ! Solar heat flux
3028      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
3029      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
3030      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
3031      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
3032      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
3033      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
3034      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
3035      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
3036
3037#if defined key_si3
3038      !                                                      ! ------------------------- !
3039      !                                                      ! Sea surface freezing temp !
3040      !                                                      ! ------------------------- !
3041      ! needed by Met Office
3042      CALL eos_fzp(tsn(:,:,1,jp_sal), sstfrz)
3043      ztmp1(:,:) = sstfrz(:,:) + rt0
3044      IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info)
3045#endif
3046      !
3047      IF (ln_timing) CALL timing_stop('sbc_cpl_snd')
3048   END SUBROUTINE sbc_cpl_snd
3049   
3050   !!======================================================================
3051END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.