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 branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/nemo_v3_6_STABLE_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 6256

Last change on this file since 6256 was 6256, checked in by frrh, 9 years ago

Add some TEMPORARY code to facilitate testing MEDUSA coupling
interface fields with UM atmosphere.

File size: 145.9 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   !!   namsbc_cpl      : coupled formulation namlist
13   !!   sbc_cpl_init    : initialisation of the coupled exchanges
14   !!   sbc_cpl_rcv     : receive fields from the atmosphere over the ocean (ocean only)
15   !!                     receive stress from the atmosphere over the ocean (ocean-ice case)
16   !!   sbc_cpl_ice_tau : receive stress from the atmosphere over ice
17   !!   sbc_cpl_ice_flx : receive fluxes from the atmosphere over ice
18   !!   sbc_cpl_snd     : send     fields to the atmosphere
19   !!----------------------------------------------------------------------
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce         ! Surface boundary condition: ocean fields
22   USE sbc_ice         ! Surface boundary condition: ice fields
23   USE sbcapr
24   USE sbcdcy          ! surface boundary condition: diurnal cycle
25   USE phycst          ! physical constants
26#if defined key_lim3
27   USE ice             ! ice variables
28#endif
29#if defined key_lim2
30   USE par_ice_2       ! ice parameters
31   USE ice_2           ! ice variables
32#endif
33   USE cpl_oasis3      ! OASIS3 coupling
34   USE geo2ocean       !
35   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev
36   USE albedo          !
37   USE in_out_manager  ! I/O manager
38   USE iom             ! NetCDF library
39   USE lib_mpp         ! distribued memory computing library
40   USE wrk_nemo        ! work arrays
41   USE timing          ! Timing
42   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
43#if defined key_oasis3 || defined key_oasis3mct
44   USE mod_oasis                    ! OASIS3-MCT module
45#endif
46   USE eosbn2
47   USE sbcrnf   , ONLY : l_rnfcpl
48#if defined key_cpl_carbon_cycle
49   USE p4zflx, ONLY : oce_co2
50#endif
51#if defined key_lim3
52   USE limthd_dh       ! for CALL lim_thd_snwblow
53#endif
54#if defined key_medusa
55! RSRH: Temporarily commented out until we get MEDUSA code included
56!   USE trc
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 sbc_ice_lim(_2).F90
64   PUBLIC   sbc_cpl_snd        ! routine called by step.F90
65   PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
66   PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).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_ts_ice = 43            ! skin temperature of sea-ice (used for melt-ponds)
112   INTEGER, PARAMETER ::   jpr_grnm   = 44            ! Greenland ice mass
113   INTEGER, PARAMETER ::   jpr_antm   = 45            ! Antarctic ice mass
114   INTEGER, PARAMETER ::   jpr_atm_pco2 = 46          ! Incoming atm CO2 flux
115   INTEGER, PARAMETER ::   jpr_atm_dust = 47          ! Incoming atm aggregate dust
116   INTEGER, PARAMETER ::   jprcv      = 47            ! total number of fields received
117
118   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
119   INTEGER, PARAMETER ::   jps_toce   =  2            ! ocean temperature
120   INTEGER, PARAMETER ::   jps_tice   =  3            ! ice   temperature
121   INTEGER, PARAMETER ::   jps_tmix   =  4            ! mixed temperature (ocean+ice)
122   INTEGER, PARAMETER ::   jps_albice =  5            ! ice   albedo
123   INTEGER, PARAMETER ::   jps_albmix =  6            ! mixed albedo
124   INTEGER, PARAMETER ::   jps_hice   =  7            ! ice  thickness
125   INTEGER, PARAMETER ::   jps_hsnw   =  8            ! snow thickness
126   INTEGER, PARAMETER ::   jps_ocx1   =  9            ! ocean current on grid 1
127   INTEGER, PARAMETER ::   jps_ocy1   = 10            !
128   INTEGER, PARAMETER ::   jps_ocz1   = 11            !
129   INTEGER, PARAMETER ::   jps_ivx1   = 12            ! ice   current on grid 1
130   INTEGER, PARAMETER ::   jps_ivy1   = 13            !
131   INTEGER, PARAMETER ::   jps_ivz1   = 14            !
132   INTEGER, PARAMETER ::   jps_co2    = 15
133   INTEGER, PARAMETER ::   jps_soce   = 16            ! ocean salinity
134   INTEGER, PARAMETER ::   jps_ssh    = 17            ! sea surface height
135   INTEGER, PARAMETER ::   jps_qsroce = 18            ! Qsr above the ocean
136   INTEGER, PARAMETER ::   jps_qnsoce = 19            ! Qns above the ocean
137   INTEGER, PARAMETER ::   jps_oemp   = 20            ! ocean freshwater budget (evap - precip)
138   INTEGER, PARAMETER ::   jps_sflx   = 21            ! salt flux
139   INTEGER, PARAMETER ::   jps_otx1   = 22            ! 2 atmosphere-ocean stress components on grid 1
140   INTEGER, PARAMETER ::   jps_oty1   = 23            !
141   INTEGER, PARAMETER ::   jps_rnf    = 24            ! runoffs
142   INTEGER, PARAMETER ::   jps_taum   = 25            ! wind stress module
143   INTEGER, PARAMETER ::   jps_fice2  = 26            ! ice fraction sent to OPA (by SAS when doing SAS-OPA coupling)
144   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
145   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
146   INTEGER, PARAMETER ::   jps_a_p    = 29            ! meltpond fraction 
147   INTEGER, PARAMETER ::   jps_ht_p   = 30            ! meltpond depth (m)
148   INTEGER, PARAMETER ::   jps_kice   = 31            ! ice surface layer thermal conductivity
149   INTEGER, PARAMETER ::   jps_sstfrz = 32            ! sea-surface freezing temperature
150   INTEGER, PARAMETER ::   jps_fice1  = 33            ! first-order ice concentration (for time-travelling ice coupling)
151   INTEGER, PARAMETER ::   jps_bio_co2 = 34           ! MEDUSA air-sea CO2 flux in
152   INTEGER, PARAMETER ::   jps_bio_dms = 35           ! MEDUSA DMS surface concentration in
153   INTEGER, PARAMETER ::   jpsnd      = 35            ! total number of fields sent
154
155   REAL(wp), PARAMETER :: dms_unit_conv = 1.0e+6      ! Coversion factor to get outgong DMS in standard units for coupling
156                                                 ! i.e. specifically nmol/L (= umol/m3)
157
158   !                                                         !!** namelist namsbc_cpl **
159   TYPE ::   FLD_C
160      CHARACTER(len = 32) ::   cldes                  ! desciption of the coupling strategy
161      CHARACTER(len = 32) ::   clcat                  ! multiple ice categories strategy
162      CHARACTER(len = 32) ::   clvref                 ! reference of vector ('spherical' or 'cartesian')
163      CHARACTER(len = 32) ::   clvor                  ! orientation of vector fields ('eastward-northward' or 'local grid')
164      CHARACTER(len = 32) ::   clvgrd                 ! grids on which is located the vector fields
165   END TYPE FLD_C
166   ! Send to the atmosphere                           !
167   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2, sn_snd_cond, sn_snd_mpnd, sn_snd_sstfrz, sn_snd_thick1
168   TYPE(FLD_C) ::   sn_snd_bio_co2, sn_snd_bio_dms                       
169
170   ! Received from the atmosphere                     !
171   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf
172   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm
173   TYPE(FLD_C) ::   sn_rcv_atm_pco2, sn_rcv_atm_dust                         
174
175   ! Other namelist parameters                        !
176   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
177   LOGICAL     ::   ln_usecplmask          !  use a coupling mask file to merge data received from several models
178                                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
179   TYPE ::   DYNARR     
180      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z3   
181   END TYPE DYNARR
182
183   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                      ! all fields recieved from the atmosphere
184
185   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix     ! ocean albedo sent to atmosphere (mix clear/overcast sky)
186
187   INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
188
189   !! Substitution
190#  include "domzgr_substitute.h90"
191#  include "vectopt_loop_substitute.h90"
192   !!----------------------------------------------------------------------
193   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
194   !! $Id$
195   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
196   !!----------------------------------------------------------------------
197
198CONTAINS
199 
200   INTEGER FUNCTION sbc_cpl_alloc()
201      !!----------------------------------------------------------------------
202      !!             ***  FUNCTION sbc_cpl_alloc  ***
203      !!----------------------------------------------------------------------
204      INTEGER :: ierr(3)
205      !!----------------------------------------------------------------------
206      ierr(:) = 0
207      !
208      ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) )
209     
210#if ! defined key_lim3 && ! defined key_lim2 && ! defined key_cice
211      ALLOCATE( a_i(jpi,jpj,1) , STAT=ierr(2) )  ! used in sbcice_if.F90 (done here as there is no sbc_ice_if_init)
212#endif
213      !ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) )
214      ! Hardwire only two models as nn_cplmodel has not been read in
215      ! from the namelist yet.
216      ALLOCATE( xcplmask(jpi,jpj,0:2) , STAT=ierr(3) )   
217      !
218      sbc_cpl_alloc = MAXVAL( ierr )
219      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
220      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
221      !
222   END FUNCTION sbc_cpl_alloc
223
224
225   SUBROUTINE sbc_cpl_init( k_ice )     
226      !!----------------------------------------------------------------------
227      !!             ***  ROUTINE sbc_cpl_init  ***
228      !!
229      !! ** Purpose :   Initialisation of send and received information from
230      !!                the atmospheric component
231      !!
232      !! ** Method  : * Read namsbc_cpl namelist
233      !!              * define the receive interface
234      !!              * define the send    interface
235      !!              * initialise the OASIS coupler
236      !!----------------------------------------------------------------------
237      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3)
238      !!
239      INTEGER ::   jn   ! dummy loop index
240      INTEGER ::   ios  ! Local integer output status for namelist read
241      INTEGER ::   inum 
242      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
243      !!
244      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick , sn_snd_crt   , sn_snd_co2,     &
245         &                  sn_snd_cond, sn_snd_mpnd  , sn_snd_sstfrz, sn_snd_thick1,                 &
246         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     &
247         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
248         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
249         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, &
250         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance,             &
251         &                  sn_snd_bio_co2, sn_snd_bio_dms,                                           &
252         &                  sn_rcv_atm_pco2, sn_rcv_atm_dust
253      !!---------------------------------------------------------------------
254      !
255      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_init')
256      !
257      CALL wrk_alloc( jpi,jpj, zacs, zaos )
258
259      ! ================================ !
260      !      Namelist informations       !
261      ! ================================ !
262
263      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
264      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
265901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
266
267      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
268      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
269902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
270      IF(lwm) WRITE ( numond, namsbc_cpl )
271
272      IF(lwp) THEN                        ! control print
273         WRITE(numout,*)
274         WRITE(numout,*)'sbc_cpl_init : namsbc_cpl namelist '
275         WRITE(numout,*)'~~~~~~~~~~~~'
276      ENDIF
277      IF( lwp .AND. ln_cpl ) THEN                        ! control print
278         WRITE(numout,*)'  received fields (mutiple ice categories)'
279         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
280         WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
281         WRITE(numout,*)'      surface stress                  = ', TRIM(sn_rcv_tau%cldes   ), ' (', TRIM(sn_rcv_tau%clcat   ), ')'
282         WRITE(numout,*)'                     - referential    = ', sn_rcv_tau%clvref
283         WRITE(numout,*)'                     - orientation    = ', sn_rcv_tau%clvor
284         WRITE(numout,*)'                     - mesh           = ', sn_rcv_tau%clvgrd
285         WRITE(numout,*)'      non-solar heat flux sensitivity = ', TRIM(sn_rcv_dqnsdt%cldes), ' (', TRIM(sn_rcv_dqnsdt%clcat), ')'
286         WRITE(numout,*)'      solar heat flux                 = ', TRIM(sn_rcv_qsr%cldes   ), ' (', TRIM(sn_rcv_qsr%clcat   ), ')'
287         WRITE(numout,*)'      non-solar heat flux             = ', TRIM(sn_rcv_qns%cldes   ), ' (', TRIM(sn_rcv_qns%clcat   ), ')'
288         WRITE(numout,*)'      freshwater budget               = ', TRIM(sn_rcv_emp%cldes   ), ' (', TRIM(sn_rcv_emp%clcat   ), ')'
289         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')'
290         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')'
291         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')'
292         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')'
293         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
294         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
295         WRITE(numout,*)'      atm pco2                        = ', TRIM(sn_rcv_atm_pco2%cldes), ' (', TRIM(sn_rcv_atm_pco2%clcat), ')'
296         WRITE(numout,*)'      atm dust                        = ', TRIM(sn_rcv_atm_dust%cldes), ' (', TRIM(sn_rcv_atm_dust%clcat), ')'
297         WRITE(numout,*)'  sent fields (multiple ice categories)'
298         WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
299         WRITE(numout,*)'      albedo                          = ', TRIM(sn_snd_alb%cldes   ), ' (', TRIM(sn_snd_alb%clcat   ), ')'
300         WRITE(numout,*)'      ice/snow thickness              = ', TRIM(sn_snd_thick%cldes ), ' (', TRIM(sn_snd_thick%clcat ), ')'
301         WRITE(numout,*)'      surface current                 = ', TRIM(sn_snd_crt%cldes   ), ' (', TRIM(sn_snd_crt%clcat   ), ')'
302         WRITE(numout,*)'                      - referential   = ', sn_snd_crt%clvref 
303         WRITE(numout,*)'                      - orientation   = ', sn_snd_crt%clvor
304         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
305         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
306         WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes   ), ' (', TRIM(sn_snd_cond%clcat   ), ')'
307         WRITE(numout,*)'      meltponds fraction & depth      = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat   ), ')'
308         WRITE(numout,*)'      sea surface freezing temp       = ', TRIM(sn_snd_sstfrz%cldes   ), ' (', TRIM(sn_snd_sstfrz%clcat   ), ')'
309
310         WRITE(numout,*)'      bio co2 flux                    = ', TRIM(sn_snd_bio_co2%cldes), ' (', TRIM(sn_snd_bio_co2%clcat), ')'
311         WRITE(numout,*)'      bio dms flux                    = ', TRIM(sn_snd_bio_dms%cldes), ' (', TRIM(sn_snd_bio_dms%clcat), ')'
312         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
313         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
314         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes
315         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction
316         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction
317         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance
318      ENDIF
319
320      !                                   ! allocate sbccpl arrays
321      !IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
322     
323      ! ================================ !
324      !   Define the receive interface   !
325      ! ================================ !
326      nrcvinfo(:) = OASIS_idle   ! needed by nrcvinfo(jpr_otx1) if we do not receive ocean stress
327
328      ! for each field: define the OASIS name                              (srcv(:)%clname)
329      !                 define receive or not from the namelist parameters (srcv(:)%laction)
330      !                 define the north fold type of lbc                  (srcv(:)%nsgn)
331
332      ! default definitions of srcv
333      srcv(:)%laction = .FALSE.   ;   srcv(:)%clgrid = 'T'   ;   srcv(:)%nsgn = 1.   ;   srcv(:)%nct = 1
334
335      !                                                      ! ------------------------- !
336      !                                                      ! ice and ocean wind stress !   
337      !                                                      ! ------------------------- !
338      !                                                           ! Name
339      srcv(jpr_otx1)%clname = 'O_OTaux1'      ! 1st ocean component on grid ONE (T or U)
340      srcv(jpr_oty1)%clname = 'O_OTauy1'      ! 2nd   -      -         -     -
341      srcv(jpr_otz1)%clname = 'O_OTauz1'      ! 3rd   -      -         -     -
342      srcv(jpr_otx2)%clname = 'O_OTaux2'      ! 1st ocean component on grid TWO (V)
343      srcv(jpr_oty2)%clname = 'O_OTauy2'      ! 2nd   -      -         -     -
344      srcv(jpr_otz2)%clname = 'O_OTauz2'      ! 3rd   -      -         -     -
345      !
346      srcv(jpr_itx1)%clname = 'O_ITaux1'      ! 1st  ice  component on grid ONE (T, F, I or U)
347      srcv(jpr_ity1)%clname = 'O_ITauy1'      ! 2nd   -      -         -     -
348      srcv(jpr_itz1)%clname = 'O_ITauz1'      ! 3rd   -      -         -     -
349      srcv(jpr_itx2)%clname = 'O_ITaux2'      ! 1st  ice  component on grid TWO (V)
350      srcv(jpr_ity2)%clname = 'O_ITauy2'      ! 2nd   -      -         -     -
351      srcv(jpr_itz2)%clname = 'O_ITauz2'      ! 3rd   -      -         -     -
352      !
353      ! Vectors: change of sign at north fold ONLY if on the local grid
354      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
355     
356      !                                                           ! Set grid and action
357      SELECT CASE( TRIM( sn_rcv_tau%clvgrd ) )      !  'T', 'U,V', 'U,V,I', 'U,V,F', 'T,I', 'T,F', or 'T,U,V'
358      CASE( 'T' ) 
359         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
360         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
361         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
362      CASE( 'U,V' ) 
363         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
364         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
365         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
366         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
367         srcv(jpr_otx1:jpr_itz2)%laction = .TRUE.     ! receive oce and ice components on both grid 1 & 2
368      CASE( 'U,V,T' )
369         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
370         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
371         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'T'        ! ice components given at T-point
372         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
373         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
374      CASE( 'U,V,I' )
375         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
376         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
377         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
378         srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
379         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
380      CASE( 'U,V,F' )
381         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'U'        ! oce components given at U-point
382         srcv(jpr_otx2:jpr_otz2)%clgrid  = 'V'        !           and           V-point
383         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
384         !srcv(jpr_otx1:jpr_otz2)%laction = .TRUE.     ! receive oce components on grid 1 & 2
385! Currently needed for HadGEM3 - but shouldn't affect anyone else for the moment
386         srcv(jpr_otx1)%laction = .TRUE. 
387         srcv(jpr_oty1)%laction = .TRUE.
388!
389         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1 only
390      CASE( 'T,I' ) 
391         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
392         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'I'        ! ice components given at I-point
393         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
394         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
395      CASE( 'T,F' ) 
396         srcv(jpr_otx1:jpr_itz2)%clgrid  = 'T'        ! oce and ice components given at T-point
397         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'F'        ! ice components given at F-point
398         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1
399         srcv(jpr_itx1:jpr_itz1)%laction = .TRUE.     ! receive ice components on grid 1
400      CASE( 'T,U,V' )
401         srcv(jpr_otx1:jpr_otz1)%clgrid  = 'T'        ! oce components given at T-point
402         srcv(jpr_itx1:jpr_itz1)%clgrid  = 'U'        ! ice components given at U-point
403         srcv(jpr_itx2:jpr_itz2)%clgrid  = 'V'        !           and           V-point
404         srcv(jpr_otx1:jpr_otz1)%laction = .TRUE.     ! receive oce components on grid 1 only
405         srcv(jpr_itx1:jpr_itz2)%laction = .TRUE.     ! receive ice components on grid 1 & 2
406      CASE default   
407         CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_tau%clvgrd' )
408      END SELECT
409      !
410      IF( TRIM( sn_rcv_tau%clvref ) == 'spherical' )   &           ! spherical: 3rd component not received
411         &     srcv( (/jpr_otz1, jpr_otz2, jpr_itz1, jpr_itz2/) )%laction = .FALSE. 
412      !
413      IF( TRIM( sn_rcv_tau%clvor  ) == 'local grid' ) THEN        ! already on local grid -> no need of the second grid
414            srcv(jpr_otx2:jpr_otz2)%laction = .FALSE. 
415            srcv(jpr_itx2:jpr_itz2)%laction = .FALSE. 
416            srcv(jpr_oty1)%clgrid = srcv(jpr_oty2)%clgrid   ! not needed but cleaner...
417            srcv(jpr_ity1)%clgrid = srcv(jpr_ity2)%clgrid   ! not needed but cleaner...
418      ENDIF
419      !
420      IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
421         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
422         srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
423         srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
424      ENDIF
425       
426      !                                                      ! ------------------------- !
427      !                                                      !    freshwater budget      !   E-P
428      !                                                      ! ------------------------- !
429      ! we suppose that atmosphere modele do not make the difference between precipiration (liquide or solid)
430      ! over ice of free ocean within the same atmospheric cell.cd
431      srcv(jpr_rain)%clname = 'OTotRain'      ! Rain = liquid precipitation
432      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation
433      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation)
434      srcv(jpr_ievp)%clname = 'OIceEvp'      ! evaporation over ice = sublimation
435      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation
436      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation
437      srcv(jpr_oemp)%clname = 'OOEvaMPr'      ! ocean water budget = ocean Evap - ocean precip
438      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
439      CASE( 'none'          )       ! nothing to do
440      CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
441      CASE( 'conservative'  )
442         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
443         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
444      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
445      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
446      END SELECT
447      !Set the number of categories for coupling of sublimation
448      IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl
449      !
450      !                                                      ! ------------------------- !
451      !                                                      !     Runoffs & Calving     !   
452      !                                                      ! ------------------------- !
453      srcv(jpr_rnf   )%clname = 'O_Runoff'
454      IF( TRIM( sn_rcv_rnf%cldes ) == 'coupled' ) THEN
455         srcv(jpr_rnf)%laction = .TRUE.
456         l_rnfcpl              = .TRUE.                      ! -> no need to read runoffs in sbcrnf
457         ln_rnf                = nn_components /= jp_iam_sas ! -> force to go through sbcrnf if not sas
458         IF(lwp) WRITE(numout,*)
459         IF(lwp) WRITE(numout,*) '   runoffs received from oasis -> force ln_rnf = ', ln_rnf
460      ENDIF
461      !
462      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE.
463      srcv(jpr_grnm  )%clname = 'OGrnmass'   ;   IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' )   srcv(jpr_grnm)%laction = .TRUE.
464      srcv(jpr_antm  )%clname = 'OAntmass'   ;   IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' )   srcv(jpr_antm)%laction = .TRUE.
465
466
467      !                                                      ! ------------------------- !
468      !                                                      !    non solar radiation    !   Qns
469      !                                                      ! ------------------------- !
470      srcv(jpr_qnsoce)%clname = 'O_QnsOce'
471      srcv(jpr_qnsice)%clname = 'O_QnsIce'
472      srcv(jpr_qnsmix)%clname = 'O_QnsMix'
473      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )
474      CASE( 'none'          )       ! nothing to do
475      CASE( 'oce only'      )   ;   srcv(               jpr_qnsoce   )%laction = .TRUE.
476      CASE( 'conservative'  )   ;   srcv( (/jpr_qnsice, jpr_qnsmix/) )%laction = .TRUE.
477      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qnsice, jpr_qnsoce/) )%laction = .TRUE.
478      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qnsmix   )%laction = .TRUE. 
479      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qns%cldes' )
480      END SELECT
481      IF( TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
482         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qns%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
483      !                                                      ! ------------------------- !
484      !                                                      !    solar radiation        !   Qsr
485      !                                                      ! ------------------------- !
486      srcv(jpr_qsroce)%clname = 'O_QsrOce'
487      srcv(jpr_qsrice)%clname = 'O_QsrIce'
488      srcv(jpr_qsrmix)%clname = 'O_QsrMix'
489      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )
490      CASE( 'none'          )       ! nothing to do
491      CASE( 'oce only'      )   ;   srcv(               jpr_qsroce   )%laction = .TRUE.
492      CASE( 'conservative'  )   ;   srcv( (/jpr_qsrice, jpr_qsrmix/) )%laction = .TRUE.
493      CASE( 'oce and ice'   )   ;   srcv( (/jpr_qsrice, jpr_qsroce/) )%laction = .TRUE.
494      CASE( 'mixed oce-ice' )   ;   srcv(               jpr_qsrmix   )%laction = .TRUE. 
495      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_qsr%cldes' )
496      END SELECT
497      IF( TRIM( sn_rcv_qsr%cldes ) == 'mixed oce-ice' .AND. jpl > 1 ) &
498         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_qsr%cldes not currently allowed to be mixed oce-ice for multi-category ice' )
499      !                                                      ! ------------------------- !
500      !                                                      !   non solar sensitivity   !   d(Qns)/d(T)
501      !                                                      ! ------------------------- !
502      srcv(jpr_dqnsdt)%clname = 'O_dQnsdT'   
503      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'coupled' )   srcv(jpr_dqnsdt)%laction = .TRUE.
504      !
505      ! non solar sensitivity mandatory for LIM ice model
506      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. k_ice /= 0 .AND. k_ice /= 4 .AND. nn_components /= jp_iam_sas ) &
507         CALL ctl_stop( 'sbc_cpl_init: sn_rcv_dqnsdt%cldes must be coupled in namsbc_cpl namelist' )
508      ! non solar sensitivity mandatory for mixed oce-ice solar radiation coupling technique
509      IF( TRIM( sn_rcv_dqnsdt%cldes ) == 'none' .AND. TRIM( sn_rcv_qns%cldes ) == 'mixed oce-ice' ) &
510         CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
511      !                                                      ! ------------------------- !
512      !                                                      !      10m wind module      !   
513      !                                                      ! ------------------------- !
514      srcv(jpr_w10m)%clname = 'O_Wind10'   ;   IF( TRIM(sn_rcv_w10m%cldes  ) == 'coupled' )   srcv(jpr_w10m)%laction = .TRUE. 
515      !
516      !                                                      ! ------------------------- !
517      !                                                      !   wind stress module      !   
518      !                                                      ! ------------------------- !
519      srcv(jpr_taum)%clname = 'O_TauMod'   ;   IF( TRIM(sn_rcv_taumod%cldes) == 'coupled' )   srcv(jpr_taum)%laction = .TRUE.
520      lhftau = srcv(jpr_taum)%laction
521
522      !                                                      ! ------------------------- !
523      !                                                      !      Atmospheric CO2      !
524      !                                                      ! ------------------------- !
525      srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
526
527
528      !                                                      ! --------------------------------------- !   
529      !                                                      ! Incoming CO2 and DUST fluxes for MEDUSA !
530      !                                                      ! --------------------------------------- ! 
531      srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
532
533      IF (TRIM(sn_rcv_atm_pco2%cldes) == 'coupled') THEN
534        srcv(jpr_atm_pco2)%laction = .TRUE.
535      END IF
536               
537      srcv(jpr_atm_dust)%clname = 'OATMDUST'   
538      IF (TRIM(sn_rcv_atm_dust%cldes) == 'coupled')  THEN
539        srcv(jpr_atm_dust)%laction = .TRUE.
540      END IF
541   
542      !                                                      ! ------------------------- !
543      !                                                      !   topmelt and botmelt     !   
544      !                                                      ! ------------------------- !
545      srcv(jpr_topm )%clname = 'OTopMlt'
546      srcv(jpr_botm )%clname = 'OBotMlt'
547      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
548         IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN
549            srcv(jpr_topm:jpr_botm)%nct = jpl
550         ELSE
551            CALL ctl_stop( 'sbc_cpl_init: sn_rcv_iceflx%clcat should always be set to yes currently' )
552         ENDIF
553         srcv(jpr_topm:jpr_botm)%laction = .TRUE.
554      ENDIF
555     
556#if defined key_cice && ! defined key_cice4
557      !                                                      ! ----------------------------- !
558      !                                                      !  sea-ice skin temperature     !   
559      !                                                      !  used in meltpond scheme      !
560      !                                                      !  May be calculated in Atm     !
561      !                                                      ! ----------------------------- !
562      srcv(jpr_ts_ice)%clname = 'OTsfIce'
563      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' ) srcv(jpr_ts_ice)%laction = .TRUE.
564      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' ) srcv(jpr_ts_ice)%nct = jpl
565      !TODO: Should there be a consistency check here?
566#endif
567
568      !                                                      ! ------------------------------- !
569      !                                                      !   OPA-SAS coupling - rcv by opa !   
570      !                                                      ! ------------------------------- !
571      srcv(jpr_sflx)%clname = 'O_SFLX'
572      srcv(jpr_fice)%clname = 'RIceFrc'
573      !
574      IF( nn_components == jp_iam_opa ) THEN    ! OPA coupled to SAS via OASIS: force received field by OPA (sent by SAS)
575         srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
576         srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
577         srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
578         srcv( (/jpr_qsroce, jpr_qnsoce, jpr_oemp, jpr_sflx, jpr_fice, jpr_otx1, jpr_oty1, jpr_taum/) )%laction = .TRUE.
579         srcv(jpr_otx1)%clgrid = 'U'        ! oce components given at U-point
580         srcv(jpr_oty1)%clgrid = 'V'        !           and           V-point
581         ! Vectors: change of sign at north fold ONLY if on the local grid
582         srcv( (/jpr_otx1,jpr_oty1/) )%nsgn = -1.
583         sn_rcv_tau%clvgrd = 'U,V'
584         sn_rcv_tau%clvor = 'local grid'
585         sn_rcv_tau%clvref = 'spherical'
586         sn_rcv_emp%cldes = 'oce only'
587         !
588         IF(lwp) THEN                        ! control print
589            WRITE(numout,*)
590            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
591            WRITE(numout,*)'               OPA component  '
592            WRITE(numout,*)
593            WRITE(numout,*)'  received fields from SAS component '
594            WRITE(numout,*)'                  ice cover '
595            WRITE(numout,*)'                  oce only EMP  '
596            WRITE(numout,*)'                  salt flux  '
597            WRITE(numout,*)'                  mixed oce-ice solar flux  '
598            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
599            WRITE(numout,*)'                  wind stress U,V on local grid and sperical coordinates '
600            WRITE(numout,*)'                  wind stress module'
601            WRITE(numout,*)
602         ENDIF
603      ENDIF
604      !                                                      ! -------------------------------- !
605      !                                                      !   OPA-SAS coupling - rcv by sas  !   
606      !                                                      ! -------------------------------- !
607      srcv(jpr_toce  )%clname = 'I_SSTSST'
608      srcv(jpr_soce  )%clname = 'I_SSSal'
609      srcv(jpr_ocx1  )%clname = 'I_OCurx1'
610      srcv(jpr_ocy1  )%clname = 'I_OCury1'
611      srcv(jpr_ssh   )%clname = 'I_SSHght'
612      srcv(jpr_e3t1st)%clname = 'I_E3T1st'   
613      srcv(jpr_fraqsr)%clname = 'I_FraQsr'   
614      !
615      IF( nn_components == jp_iam_sas ) THEN
616         IF( .NOT. ln_cpl ) srcv(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
617         IF( .NOT. ln_cpl ) srcv(:)%clgrid  = 'T'       ! force default definition in case of opa <-> sas coupling
618         IF( .NOT. ln_cpl ) srcv(:)%nsgn    = 1.        ! force default definition in case of opa <-> sas coupling
619         srcv( (/jpr_toce, jpr_soce, jpr_ssh, jpr_fraqsr, jpr_ocx1, jpr_ocy1/) )%laction = .TRUE.
620         srcv( jpr_e3t1st )%laction = lk_vvl
621         srcv(jpr_ocx1)%clgrid = 'U'        ! oce components given at U-point
622         srcv(jpr_ocy1)%clgrid = 'V'        !           and           V-point
623         ! Vectors: change of sign at north fold ONLY if on the local grid
624         srcv(jpr_ocx1:jpr_ocy1)%nsgn = -1.
625         ! Change first letter to couple with atmosphere if already coupled OPA
626         ! this is nedeed as each variable name used in the namcouple must be unique:
627         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere
628         DO jn = 1, jprcv
629            IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname))
630         END DO
631         !
632         IF(lwp) THEN                        ! control print
633            WRITE(numout,*)
634            WRITE(numout,*)'               Special conditions for SAS-OPA coupling  '
635            WRITE(numout,*)'               SAS component  '
636            WRITE(numout,*)
637            IF( .NOT. ln_cpl ) THEN
638               WRITE(numout,*)'  received fields from OPA component '
639            ELSE
640               WRITE(numout,*)'  Additional received fields from OPA component : '
641            ENDIF
642            WRITE(numout,*)'               sea surface temperature (Celcius) '
643            WRITE(numout,*)'               sea surface salinity ' 
644            WRITE(numout,*)'               surface currents ' 
645            WRITE(numout,*)'               sea surface height ' 
646            WRITE(numout,*)'               thickness of first ocean T level '       
647            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
648            WRITE(numout,*)
649         ENDIF
650      ENDIF
651     
652      ! =================================================== !
653      ! Allocate all parts of frcv used for received fields !
654      ! =================================================== !
655      DO jn = 1, jprcv
656         IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) )
657      END DO
658      ! Allocate taum part of frcv which is used even when not received as coupling field
659      IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) )
660      ! Allocate w10m part of frcv which is used even when not received as coupling field
661      IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) )
662      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field
663      IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) )
664      IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) )
665      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
666      IF( k_ice /= 0 ) THEN
667         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) )
668         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) )
669      END IF
670
671      ! ================================ !
672      !     Define the send interface    !
673      ! ================================ !
674      ! for each field: define the OASIS name                           (ssnd(:)%clname)
675      !                 define send or not from the namelist parameters (ssnd(:)%laction)
676      !                 define the north fold type of lbc               (ssnd(:)%nsgn)
677     
678      ! default definitions of nsnd
679      ssnd(:)%laction = .FALSE.   ;   ssnd(:)%clgrid = 'T'   ;   ssnd(:)%nsgn = 1.  ; ssnd(:)%nct = 1
680         
681      !                                                      ! ------------------------- !
682      !                                                      !    Surface temperature    !
683      !                                                      ! ------------------------- !
684      ssnd(jps_toce)%clname = 'O_SSTSST'
685      ssnd(jps_tice)%clname = 'OTepIce'
686      ssnd(jps_tmix)%clname = 'O_TepMix'
687      SELECT CASE( TRIM( sn_snd_temp%cldes ) )
688      CASE( 'none'                                 )       ! nothing to do
689      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE.
690      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice')
691         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE.
692         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl
693      CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE.
694      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' )
695      END SELECT
696           
697      !                                                      ! ------------------------- !
698      !                                                      !          Albedo           !
699      !                                                      ! ------------------------- !
700      ssnd(jps_albice)%clname = 'O_AlbIce' 
701      ssnd(jps_albmix)%clname = 'O_AlbMix'
702      SELECT CASE( TRIM( sn_snd_alb%cldes ) )
703      CASE( 'none'                 )     ! nothing to do
704      CASE( 'ice' , 'weighted ice' )   ; ssnd(jps_albice)%laction = .TRUE.
705      CASE( 'mixed oce-ice'        )   ; ssnd(jps_albmix)%laction = .TRUE.
706      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_alb%cldes' )
707      END SELECT
708      !
709      ! Need to calculate oceanic albedo if
710      !     1. sending mixed oce-ice albedo or
711      !     2. receiving mixed oce-ice solar radiation
712      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN
713         CALL albedo_oce( zaos, zacs )
714         ! Due to lack of information on nebulosity : mean clear/overcast sky
715         albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5
716      ENDIF
717
718      !                                                      ! ------------------------- !
719      !                                                      !  Ice fraction & Thickness
720      !                                                      ! ------------------------- !
721      ssnd(jps_fice)%clname = 'OIceFrc'
722      ssnd(jps_hice)%clname = 'OIceTck'
723      ssnd(jps_hsnw)%clname = 'OSnwTck'
724      ssnd(jps_a_p)%clname  = 'OPndFrc'
725      ssnd(jps_ht_p)%clname = 'OPndTck'
726      ssnd(jps_fice1)%clname = 'OIceFrd'
727      IF( k_ice /= 0 ) THEN
728         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case)
729         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used
730                                                     ! in producing atmos-to-ice fluxes
731! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now
732         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl
733         IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = jpl
734      ENDIF
735     
736      SELECT CASE ( TRIM( sn_snd_thick%cldes ) )
737      CASE( 'none'         )       ! nothing to do
738      CASE( 'ice and snow' ) 
739         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
740         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN
741            ssnd(jps_hice:jps_hsnw)%nct = jpl
742         ENDIF
743      CASE ( 'weighted ice and snow' ) 
744         ssnd(jps_hice:jps_hsnw)%laction = .TRUE.
745         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = jpl
746      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' )
747      END SELECT
748
749      !                                                      ! ------------------------- !
750      !                                                      ! Ice Meltponds             !
751      !                                                      ! ------------------------- !
752#if defined key_cice && ! defined key_cice4
753      ! Meltponds only CICE5
754      ssnd(jps_a_p)%clname = 'OPndFrc'   
755      ssnd(jps_ht_p)%clname = 'OPndTck'   
756      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) )
757      CASE ( 'none' )
758         ssnd(jps_a_p)%laction = .FALSE.
759         ssnd(jps_ht_p)%laction = .FALSE.
760      CASE ( 'ice only' ) 
761         ssnd(jps_a_p)%laction = .TRUE.
762         ssnd(jps_ht_p)%laction = .TRUE.
763         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
764            ssnd(jps_a_p)%nct = jpl
765            ssnd(jps_ht_p)%nct = jpl
766         ELSE
767            IF ( jpl > 1 ) THEN
768               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' )
769            ENDIF
770         ENDIF
771      CASE ( 'weighted ice' ) 
772         ssnd(jps_a_p)%laction = .TRUE.
773         ssnd(jps_ht_p)%laction = .TRUE.
774         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN
775            ssnd(jps_a_p)%nct = jpl 
776            ssnd(jps_ht_p)%nct = jpl 
777         ENDIF
778      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes' )
779      END SELECT
780#else
781      IF( TRIM( sn_snd_mpnd%cldes /= 'none' ) THEN
782         CALL ctl_stop('Meltponds can only be used with CICEv5')
783      ENDIF
784#endif
785
786      !                                                      ! ------------------------- !
787      !                                                      !      Surface current      !
788      !                                                      ! ------------------------- !
789      !        ocean currents              !            ice velocities
790      ssnd(jps_ocx1)%clname = 'O_OCurx1'   ;   ssnd(jps_ivx1)%clname = 'O_IVelx1'
791      ssnd(jps_ocy1)%clname = 'O_OCury1'   ;   ssnd(jps_ivy1)%clname = 'O_IVely1'
792      ssnd(jps_ocz1)%clname = 'O_OCurz1'   ;   ssnd(jps_ivz1)%clname = 'O_IVelz1'
793      !
794      ssnd(jps_ocx1:jps_ivz1)%nsgn = -1.   ! vectors: change of the sign at the north fold
795
796      IF( sn_snd_crt%clvgrd == 'U,V' ) THEN
797         ssnd(jps_ocx1)%clgrid = 'U' ; ssnd(jps_ocy1)%clgrid = 'V'
798      ELSE IF( sn_snd_crt%clvgrd /= 'T' ) THEN 
799         CALL ctl_stop( 'sn_snd_crt%clvgrd must be equal to T' )
800         ssnd(jps_ocx1:jps_ivz1)%clgrid  = 'T'      ! all oce and ice components on the same unique grid
801      ENDIF
802      ssnd(jps_ocx1:jps_ivz1)%laction = .TRUE.   ! default: all are send
803      IF( TRIM( sn_snd_crt%clvref ) == 'spherical' )   ssnd( (/jps_ocz1, jps_ivz1/) )%laction = .FALSE. 
804      IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) ssnd(jps_ocx1:jps_ivz1)%nsgn = 1.
805      SELECT CASE( TRIM( sn_snd_crt%cldes ) )
806      CASE( 'none'                 )   ;   ssnd(jps_ocx1:jps_ivz1)%laction = .FALSE.
807      CASE( 'oce only'             )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
808      CASE( 'weighted oce and ice' )   !   nothing to do
809      CASE( 'mixed oce-ice'        )   ;   ssnd(jps_ivx1:jps_ivz1)%laction = .FALSE.
810      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_crt%cldes' )
811      END SELECT
812
813      !                                                      ! ------------------------- !
814      !                                                      !          CO2 flux         !
815      !                                                      ! ------------------------- !
816      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE.
817      !
818
819      ! MEDUSA
820      ssnd(jps_bio_co2)%clname = 'OMED_CO2'
821      IF( TRIM(sn_snd_bio_co2%cldes) == 'coupled' )    ssnd(jps_bio_co2)%laction = .TRUE.
822     
823      !                                                      ! ------------------------- !
824      !                                                      !          DMS field        !
825      !                                                      ! ------------------------- !
826      ! MEDUSA
827      ssnd(jps_bio_dms)%clname = 'OMED_DMS'
828      IF( TRIM(sn_snd_bio_dms%cldes) == 'coupled' )    ssnd(jps_bio_dms)%laction = .TRUE.
829
830     
831      !                                                      ! ------------------------- !
832      !                                                      ! Sea surface freezing temp !
833      !                                                      ! ------------------------- !
834      ssnd(jps_sstfrz)%clname = 'O_SSTFrz' ; IF( TRIM(sn_snd_sstfrz%cldes) == 'coupled' )  ssnd(jps_sstfrz)%laction = .TRUE.
835      !
836      !                                                      ! ------------------------- !
837      !                                                      !    Ice conductivity       !
838      !                                                      ! ------------------------- !
839      ! Note that ultimately we will move to passing an ocean effective conductivity as well so there
840      ! will be some changes to the parts of the code which currently relate only to ice conductivity
841      ssnd(jps_kice )%clname = 'OIceKn'
842      SELECT CASE ( TRIM( sn_snd_cond%cldes ) )
843      CASE ( 'none' )
844         ssnd(jps_kice)%laction = .FALSE.
845      CASE ( 'ice only' )
846         ssnd(jps_kice)%laction = .TRUE.
847         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN
848            ssnd(jps_kice)%nct = jpl
849         ELSE
850            IF ( jpl > 1 ) THEN
851               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' )
852            ENDIF
853         ENDIF
854      CASE ( 'weighted ice' )
855         ssnd(jps_kice)%laction = .TRUE.
856         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = jpl
857      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes' )
858      END SELECT
859      !
860     
861
862      !                                                      ! ------------------------------- !
863      !                                                      !   OPA-SAS coupling - snd by opa !   
864      !                                                      ! ------------------------------- !
865      ssnd(jps_ssh   )%clname = 'O_SSHght' 
866      ssnd(jps_soce  )%clname = 'O_SSSal' 
867      ssnd(jps_e3t1st)%clname = 'O_E3T1st'   
868      ssnd(jps_fraqsr)%clname = 'O_FraQsr'
869      !
870      IF( nn_components == jp_iam_opa ) THEN
871         ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
872         ssnd( (/jps_toce, jps_soce, jps_ssh, jps_fraqsr, jps_ocx1, jps_ocy1/) )%laction = .TRUE.
873         ssnd( jps_e3t1st )%laction = lk_vvl
874         ! vector definition: not used but cleaner...
875         ssnd(jps_ocx1)%clgrid  = 'U'        ! oce components given at U-point
876         ssnd(jps_ocy1)%clgrid  = 'V'        !           and           V-point
877         sn_snd_crt%clvgrd = 'U,V'
878         sn_snd_crt%clvor = 'local grid'
879         sn_snd_crt%clvref = 'spherical'
880         !
881         IF(lwp) THEN                        ! control print
882            WRITE(numout,*)
883            WRITE(numout,*)'  sent fields to SAS component '
884            WRITE(numout,*)'               sea surface temperature (T before, Celcius) '
885            WRITE(numout,*)'               sea surface salinity ' 
886            WRITE(numout,*)'               surface currents U,V on local grid and spherical coordinates' 
887            WRITE(numout,*)'               sea surface height ' 
888            WRITE(numout,*)'               thickness of first ocean T level '       
889            WRITE(numout,*)'               fraction of solar net radiation absorbed in the first ocean level'
890            WRITE(numout,*)
891         ENDIF
892      ENDIF
893      !                                                      ! ------------------------------- !
894      !                                                      !   OPA-SAS coupling - snd by sas !   
895      !                                                      ! ------------------------------- !
896      ssnd(jps_sflx  )%clname = 'I_SFLX'     
897      ssnd(jps_fice2 )%clname = 'IIceFrc'
898      ssnd(jps_qsroce)%clname = 'I_QsrOce'   
899      ssnd(jps_qnsoce)%clname = 'I_QnsOce'   
900      ssnd(jps_oemp  )%clname = 'IOEvaMPr' 
901      ssnd(jps_otx1  )%clname = 'I_OTaux1'   
902      ssnd(jps_oty1  )%clname = 'I_OTauy1'   
903      ssnd(jps_rnf   )%clname = 'I_Runoff'   
904      ssnd(jps_taum  )%clname = 'I_TauMod'   
905      !
906      IF( nn_components == jp_iam_sas ) THEN
907         IF( .NOT. ln_cpl ) ssnd(:)%laction = .FALSE.   ! force default definition in case of opa <-> sas coupling
908         ssnd( (/jps_qsroce, jps_qnsoce, jps_oemp, jps_fice2, jps_sflx, jps_otx1, jps_oty1, jps_taum/) )%laction = .TRUE.
909         !
910         ! Change first letter to couple with atmosphere if already coupled with sea_ice
911         ! this is nedeed as each variable name used in the namcouple must be unique:
912         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere
913         DO jn = 1, jpsnd
914            IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname))
915         END DO
916         !
917         IF(lwp) THEN                        ! control print
918            WRITE(numout,*)
919            IF( .NOT. ln_cpl ) THEN
920               WRITE(numout,*)'  sent fields to OPA component '
921            ELSE
922               WRITE(numout,*)'  Additional sent fields to OPA component : '
923            ENDIF
924            WRITE(numout,*)'                  ice cover '
925            WRITE(numout,*)'                  oce only EMP  '
926            WRITE(numout,*)'                  salt flux  '
927            WRITE(numout,*)'                  mixed oce-ice solar flux  '
928            WRITE(numout,*)'                  mixed oce-ice non solar flux  '
929            WRITE(numout,*)'                  wind stress U,V components'
930            WRITE(numout,*)'                  wind stress module'
931         ENDIF
932      ENDIF
933
934      !
935      ! ================================ !
936      !   initialisation of the coupler  !
937      ! ================================ !
938
939      CALL cpl_define(jprcv, jpsnd, nn_cplmodel)
940     
941      IF (ln_usecplmask) THEN
942         xcplmask(:,:,:) = 0.
943         CALL iom_open( 'cplmask', inum )
944         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   &
945            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) )
946         CALL iom_close( inum )
947      ELSE
948         xcplmask(:,:,:) = 1.
949      ENDIF
950      xcplmask(:,:,0) = 1. - SUM( xcplmask(:,:,1:nn_cplmodel), dim = 3 )
951      !
952      ncpl_qsr_freq = cpl_freq( 'O_QsrOce' ) + cpl_freq( 'O_QsrMix' ) + cpl_freq( 'I_QsrOce' ) + cpl_freq( 'I_QsrMix' )
953      IF( ln_dm2dc .AND. ln_cpl .AND. ncpl_qsr_freq /= 86400 )   &
954         &   CALL ctl_stop( 'sbc_cpl_init: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' )
955      ncpl_qsr_freq = 86400 / ncpl_qsr_freq
956
957      IF( ln_coupled_iceshelf_fluxes ) THEN
958          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
959          ! more complicated could be done if required.
960          greenland_icesheet_mask = 0.0
961          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0
962          antarctica_icesheet_mask = 0.0
963          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0
964
965          ! initialise other variables
966          greenland_icesheet_mass_array(:,:) = 0.0
967          antarctica_icesheet_mass_array(:,:) = 0.0
968
969          IF( .not. ln_rstart ) THEN
970             greenland_icesheet_mass = 0.0 
971             greenland_icesheet_mass_rate_of_change = 0.0 
972             greenland_icesheet_timelapsed = 0.0
973             antarctica_icesheet_mass = 0.0 
974             antarctica_icesheet_mass_rate_of_change = 0.0 
975             antarctica_icesheet_timelapsed = 0.0
976          ENDIF
977
978      ENDIF
979
980      CALL wrk_dealloc( jpi,jpj, zacs, zaos )
981      !
982      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_init')
983      !
984   END SUBROUTINE sbc_cpl_init
985
986
987   SUBROUTINE sbc_cpl_rcv( kt, k_fsbc, k_ice )     
988      !!----------------------------------------------------------------------
989      !!             ***  ROUTINE sbc_cpl_rcv  ***
990      !!
991      !! ** Purpose :   provide the stress over the ocean and, if no sea-ice,
992      !!                provide the ocean heat and freshwater fluxes.
993      !!
994      !! ** Method  : - Receive all the atmospheric fields (stored in frcv array). called at each time step.
995      !!                OASIS controls if there is something do receive or not. nrcvinfo contains the info
996      !!                to know if the field was really received or not
997      !!
998      !!              --> If ocean stress was really received:
999      !!
1000      !!                  - transform the received ocean stress vector from the received
1001      !!                 referential and grid into an atmosphere-ocean stress in
1002      !!                 the (i,j) ocean referencial and at the ocean velocity point.
1003      !!                    The received stress are :
1004      !!                     - defined by 3 components (if cartesian coordinate)
1005      !!                            or by 2 components (if spherical)
1006      !!                     - oriented along geographical   coordinate (if eastward-northward)
1007      !!                            or  along the local grid coordinate (if local grid)
1008      !!                     - given at U- and V-point, resp.   if received on 2 grids
1009      !!                            or at T-point               if received on 1 grid
1010      !!                    Therefore and if necessary, they are successively
1011      !!                  processed in order to obtain them
1012      !!                     first  as  2 components on the sphere
1013      !!                     second as  2 components oriented along the local grid
1014      !!                     third  as  2 components on the U,V grid
1015      !!
1016      !!              -->
1017      !!
1018      !!              - In 'ocean only' case, non solar and solar ocean heat fluxes
1019      !!             and total ocean freshwater fluxes 
1020      !!
1021      !! ** Method  :   receive all fields from the atmosphere and transform
1022      !!              them into ocean surface boundary condition fields
1023      !!
1024      !! ** Action  :   update  utau, vtau   ocean stress at U,V grid
1025      !!                        taum         wind stress module at T-point
1026      !!                        wndm         wind speed  module at T-point over free ocean or leads in presence of sea-ice
1027      !!                        qns          non solar heat fluxes including emp heat content    (ocean only case)
1028      !!                                     and the latent heat flux of solid precip. melting
1029      !!                        qsr          solar ocean heat fluxes   (ocean only case)
1030      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
1031      !!----------------------------------------------------------------------
1032      INTEGER, INTENT(in)           ::   kt          ! ocean model time step index
1033      INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation
1034      INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3)
1035
1036      !!
1037      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module??
1038      INTEGER  ::   ji, jj, jl, jn         ! dummy loop indices
1039      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
1040      INTEGER  ::   ikchoix
1041      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
1042      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in
1043      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b
1044      REAL(wp) ::   zmask_sum, zepsilon     
1045      REAL(wp) ::   zcoef                  ! temporary scalar
1046      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
1047      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient
1048      REAL(wp) ::   zzx, zzy               ! temporary variables
1049      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2
1050      !!----------------------------------------------------------------------
1051
1052      ! RSRH temporary arrays for testing, just to recieve incoming MEDUSA related fields
1053      ! until we know where they need to go.
1054      REAL(wp), ALLOCATABLE :: atm_pco2(:,:)
1055      REAL(wp), ALLOCATABLE :: atm_dust(:,:)
1056
1057      !
1058      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
1059      !
1060      CALL wrk_alloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2 )
1061      !
1062      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1063      !
1064      !                                                      ! ======================================================= !
1065      !                                                      ! Receive all the atmos. fields (including ice information)
1066      !                                                      ! ======================================================= !
1067      isec = ( kt - nit000 ) * NINT( rdttra(1) )                ! date of exchanges
1068      DO jn = 1, jprcv                                          ! received fields sent by the atmosphere
1069         IF( srcv(jn)%laction )   CALL cpl_rcv( jn, isec, frcv(jn)%z3, xcplmask(:,:,1:nn_cplmodel), nrcvinfo(jn) )
1070      END DO
1071
1072      !                                                      ! ========================= !
1073      IF( srcv(jpr_otx1)%laction ) THEN                      !  ocean stress components  !
1074         !                                                   ! ========================= !
1075         ! define frcv(jpr_otx1)%z3(:,:,1) and frcv(jpr_oty1)%z3(:,:,1): stress at U/V point along model grid
1076         ! => need to be done only when we receive the field
1077         IF(  nrcvinfo(jpr_otx1) == OASIS_Rcv ) THEN
1078            !
1079            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1080               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1081               !
1082               CALL geo2oce( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), frcv(jpr_otz1)%z3(:,:,1),   &
1083                  &          srcv(jpr_otx1)%clgrid, ztx, zty )
1084               frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1085               frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1086               !
1087               IF( srcv(jpr_otx2)%laction ) THEN
1088                  CALL geo2oce( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), frcv(jpr_otz2)%z3(:,:,1),   &
1089                     &          srcv(jpr_otx2)%clgrid, ztx, zty )
1090                  frcv(jpr_otx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1091                  frcv(jpr_oty2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1092               ENDIF
1093               !
1094            ENDIF
1095            !
1096            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1097               !                                                       ! (geographical to local grid -> rotate the components)
1098               IF( srcv(jpr_otx1)%clgrid == 'U' .AND. (.NOT. srcv(jpr_otx2)%laction) ) THEN
1099                  ! Temporary code for HadGEM3 - will be removed eventually.
1100        ! Only applies when we have only taux on U grid and tauy on V grid
1101             DO jj=2,jpjm1
1102                DO ji=2,jpim1
1103                     ztx(ji,jj)=0.25*vmask(ji,jj,1)                &
1104                        *(frcv(jpr_otx1)%z3(ji,jj,1)+frcv(jpr_otx1)%z3(ji-1,jj,1)    &
1105                        +frcv(jpr_otx1)%z3(ji,jj+1,1)+frcv(jpr_otx1)%z3(ji-1,jj+1,1))
1106                     zty(ji,jj)=0.25*umask(ji,jj,1)                &
1107                        *(frcv(jpr_oty1)%z3(ji,jj,1)+frcv(jpr_oty1)%z3(ji+1,jj,1)    &
1108                        +frcv(jpr_oty1)%z3(ji,jj-1,1)+frcv(jpr_oty1)%z3(ji+1,jj-1,1))
1109                ENDDO
1110             ENDDO
1111                   
1112             ikchoix = 1
1113             CALL repcmo (frcv(jpr_otx1)%z3(:,:,1),zty,ztx,frcv(jpr_oty1)%z3(:,:,1),ztx2,zty2,ikchoix)
1114             CALL lbc_lnk (ztx2,'U', -1. )
1115             CALL lbc_lnk (zty2,'V', -1. )
1116             frcv(jpr_otx1)%z3(:,:,1)=ztx2(:,:)
1117             frcv(jpr_oty1)%z3(:,:,1)=zty2(:,:)
1118          ELSE
1119             CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->i', ztx )   
1120             frcv(jpr_otx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1121             IF( srcv(jpr_otx2)%laction ) THEN
1122                CALL rot_rep( frcv(jpr_otx2)%z3(:,:,1), frcv(jpr_oty2)%z3(:,:,1), srcv(jpr_otx2)%clgrid, 'en->j', zty )   
1123             ELSE
1124                CALL rot_rep( frcv(jpr_otx1)%z3(:,:,1), frcv(jpr_oty1)%z3(:,:,1), srcv(jpr_otx1)%clgrid, 'en->j', zty ) 
1125             ENDIF
1126          frcv(jpr_oty1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 2nd grid 
1127               ENDIF
1128            ENDIF
1129            !                             
1130            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN
1131               DO jj = 2, jpjm1                                          ! T ==> (U,V)
1132                  DO ji = fs_2, fs_jpim1   ! vector opt.
1133                     frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) )
1134                     frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) )
1135                  END DO
1136               END DO
1137               CALL lbc_lnk( frcv(jpr_otx1)%z3(:,:,1), 'U',  -1. )   ;   CALL lbc_lnk( frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. )
1138            ENDIF
1139            llnewtx = .TRUE.
1140         ELSE
1141            llnewtx = .FALSE.
1142         ENDIF
1143         !                                                   ! ========================= !
1144      ELSE                                                   !   No dynamical coupling   !
1145         !                                                   ! ========================= !
1146         frcv(jpr_otx1)%z3(:,:,1) = 0.e0                               ! here simply set to zero
1147         frcv(jpr_oty1)%z3(:,:,1) = 0.e0                               ! an external read in a file can be added instead
1148         llnewtx = .TRUE.
1149         !
1150      ENDIF
1151      !                                                      ! ========================= !
1152      !                                                      !    wind stress module     !   (taum)
1153      !                                                      ! ========================= !
1154      !
1155      IF( .NOT. srcv(jpr_taum)%laction ) THEN                    ! compute wind stress module from its components if not received
1156         ! => need to be done only when otx1 was changed
1157         IF( llnewtx ) THEN
1158!CDIR NOVERRCHK
1159            DO jj = 2, jpjm1
1160!CDIR NOVERRCHK
1161               DO ji = fs_2, fs_jpim1   ! vect. opt.
1162                  zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1)
1163                  zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1)
1164                  frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy )
1165               END DO
1166            END DO
1167            CALL lbc_lnk( frcv(jpr_taum)%z3(:,:,1), 'T', 1. )
1168            llnewtau = .TRUE.
1169         ELSE
1170            llnewtau = .FALSE.
1171         ENDIF
1172      ELSE
1173         llnewtau = nrcvinfo(jpr_taum) == OASIS_Rcv
1174         ! Stress module can be negative when received (interpolation problem)
1175         IF( llnewtau ) THEN
1176            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
1177         ENDIF
1178      ENDIF
1179      !
1180      !                                                      ! ========================= !
1181      !                                                      !      10 m wind speed      !   (wndm)
1182      !                                                      ! ========================= !
1183      !
1184      IF( .NOT. srcv(jpr_w10m)%laction ) THEN                    ! compute wind spreed from wind stress module if not received 
1185         ! => need to be done only when taumod was changed
1186         IF( llnewtau ) THEN
1187            zcoef = 1. / ( zrhoa * zcdrag ) 
1188!CDIR NOVERRCHK
1189            DO jj = 1, jpj
1190!CDIR NOVERRCHK
1191               DO ji = 1, jpi 
1192                  frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef )
1193               END DO
1194            END DO
1195         ENDIF
1196      ENDIF
1197
1198      ! u(v)tau and taum will be modified by ice model
1199      ! -> need to be reset before each call of the ice/fsbc     
1200      IF( MOD( kt-1, k_fsbc ) == 0 ) THEN
1201         !
1202         IF( ln_mixcpl ) THEN
1203            utau(:,:) = utau(:,:) * xcplmask(:,:,0) + frcv(jpr_otx1)%z3(:,:,1) * zmsk(:,:)
1204            vtau(:,:) = vtau(:,:) * xcplmask(:,:,0) + frcv(jpr_oty1)%z3(:,:,1) * zmsk(:,:)
1205            taum(:,:) = taum(:,:) * xcplmask(:,:,0) + frcv(jpr_taum)%z3(:,:,1) * zmsk(:,:)
1206            wndm(:,:) = wndm(:,:) * xcplmask(:,:,0) + frcv(jpr_w10m)%z3(:,:,1) * zmsk(:,:)
1207         ELSE
1208            utau(:,:) = frcv(jpr_otx1)%z3(:,:,1)
1209            vtau(:,:) = frcv(jpr_oty1)%z3(:,:,1)
1210            taum(:,:) = frcv(jpr_taum)%z3(:,:,1)
1211            wndm(:,:) = frcv(jpr_w10m)%z3(:,:,1)
1212         ENDIF
1213         CALL iom_put( "taum_oce", taum )   ! output wind stress module
1214         
1215      ENDIF
1216
1217#if defined key_cpl_carbon_cycle
1218      !                                                      ! ================== !
1219      !                                                      ! atmosph. CO2 (ppm) !
1220      !                                                      ! ================== !
1221      IF( srcv(jpr_co2)%laction )   atm_co2(:,:) = frcv(jpr_co2)%z3(:,:,1)
1222#endif
1223
1224#if defined key_cice && ! defined key_cice4
1225      !  ! Sea ice surface skin temp:
1226      IF( srcv(jpr_ts_ice)%laction ) THEN
1227        DO jl = 1, jpl
1228          DO jj = 1, jpj
1229            DO ji = 1, jpi
1230              IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) > 0.0) THEN
1231                tsfc_ice(ji,jj,jl) = 0.0
1232              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) < -60.0) THEN
1233                tsfc_ice(ji,jj,jl) = -60.0
1234              ELSE
1235                tsfc_ice(ji,jj,jl) = frcv(jpr_ts_ice)%z3(ji,jj,jl)
1236              ENDIF
1237            END DO
1238          END DO
1239        END DO
1240      ENDIF
1241#endif
1242
1243#if defined key_medusa
1244      ! RSRH Allocate temporary arrays to receive incoming fields during testing
1245      ALLOCATE(atm_pco2(jpi,jpj))
1246      ALLOCATE(atm_dust(jpi,jpj))
1247
1248      IF( srcv(jpr_atm_pco2)%laction) atm_pco2(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
1249      IF( srcv(jpr_atm_dust)%laction) atm_dust(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
1250     
1251      ! RSRH Deallocate temporary arrays.
1252      DEALLOCATE(atm_pco2)
1253      DEALLOCATE(atm_dust)
1254#endif
1255
1256
1257
1258
1259      !  Fields received by SAS when OASIS coupling
1260      !  (arrays no more filled at sbcssm stage)
1261      !                                                      ! ================== !
1262      !                                                      !        SSS         !
1263      !                                                      ! ================== !
1264      IF( srcv(jpr_soce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1265         sss_m(:,:) = frcv(jpr_soce)%z3(:,:,1)
1266         CALL iom_put( 'sss_m', sss_m )
1267      ENDIF
1268      !                                               
1269      !                                                      ! ================== !
1270      !                                                      !        SST         !
1271      !                                                      ! ================== !
1272      IF( srcv(jpr_toce)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1273         sst_m(:,:) = frcv(jpr_toce)%z3(:,:,1)
1274         IF( srcv(jpr_soce)%laction .AND. ln_useCT ) THEN    ! make sure that sst_m is the potential temperature
1275            sst_m(:,:) = eos_pt_from_ct( sst_m(:,:), sss_m(:,:) )
1276         ENDIF
1277      ENDIF
1278      !                                                      ! ================== !
1279      !                                                      !        SSH         !
1280      !                                                      ! ================== !
1281      IF( srcv(jpr_ssh )%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1282         ssh_m(:,:) = frcv(jpr_ssh )%z3(:,:,1)
1283         CALL iom_put( 'ssh_m', ssh_m )
1284      ENDIF
1285      !                                                      ! ================== !
1286      !                                                      !  surface currents  !
1287      !                                                      ! ================== !
1288      IF( srcv(jpr_ocx1)%laction ) THEN                      ! received by sas in case of opa <-> sas coupling
1289         ssu_m(:,:) = frcv(jpr_ocx1)%z3(:,:,1)
1290         ub (:,:,1) = ssu_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1291         CALL iom_put( 'ssu_m', ssu_m )
1292      ENDIF
1293      IF( srcv(jpr_ocy1)%laction ) THEN
1294         ssv_m(:,:) = frcv(jpr_ocy1)%z3(:,:,1)
1295         vb (:,:,1) = ssv_m(:,:)                             ! will be used in sbcice_lim in the call of lim_sbc_tau
1296         CALL iom_put( 'ssv_m', ssv_m )
1297      ENDIF
1298      !                                                      ! ======================== !
1299      !                                                      !  first T level thickness !
1300      !                                                      ! ======================== !
1301      IF( srcv(jpr_e3t1st )%laction ) THEN                   ! received by sas in case of opa <-> sas coupling
1302         e3t_m(:,:) = frcv(jpr_e3t1st )%z3(:,:,1)
1303         CALL iom_put( 'e3t_m', e3t_m(:,:) )
1304      ENDIF
1305      !                                                      ! ================================ !
1306      !                                                      !  fraction of solar net radiation !
1307      !                                                      ! ================================ !
1308      IF( srcv(jpr_fraqsr)%laction ) THEN                    ! received by sas in case of opa <-> sas coupling
1309         frq_m(:,:) = frcv(jpr_fraqsr)%z3(:,:,1)
1310         CALL iom_put( 'frq_m', frq_m )
1311      ENDIF
1312     
1313      !                                                      ! ========================= !
1314      IF( k_ice <= 1 .AND. MOD( kt-1, k_fsbc ) == 0 ) THEN   !  heat & freshwater fluxes ! (Ocean only case)
1315         !                                                   ! ========================= !
1316         !
1317         !                                                       ! total freshwater fluxes over the ocean (emp)
1318         IF( srcv(jpr_oemp)%laction .OR. srcv(jpr_rain)%laction ) THEN
1319            SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
1320            CASE( 'conservative' )
1321               zemp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ( frcv(jpr_rain)%z3(:,:,1) + frcv(jpr_snow)%z3(:,:,1) )
1322            CASE( 'oce only', 'oce and ice' )
1323               zemp(:,:) = frcv(jpr_oemp)%z3(:,:,1)
1324            CASE default
1325               CALL ctl_stop( 'sbc_cpl_rcv: wrong definition of sn_rcv_emp%cldes' )
1326            END SELECT
1327         ELSE
1328            zemp(:,:) = 0._wp
1329         ENDIF
1330         !
1331         !                                                        ! runoffs and calving (added in emp)
1332         IF( srcv(jpr_rnf)%laction )     rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1333         IF( srcv(jpr_cal)%laction )     zemp(:,:) = zemp(:,:) - frcv(jpr_cal)%z3(:,:,1)
1334         
1335         IF( ln_mixcpl ) THEN   ;   emp(:,:) = emp(:,:) * xcplmask(:,:,0) + zemp(:,:) * zmsk(:,:)
1336         ELSE                   ;   emp(:,:) =                              zemp(:,:)
1337         ENDIF
1338         !
1339         !                                                       ! non solar heat flux over the ocean (qns)
1340         IF(      srcv(jpr_qnsoce)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsoce)%z3(:,:,1)
1341         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
1342         ELSE                                       ;   zqns(:,:) = 0._wp
1343         END IF
1344         ! update qns over the free ocean with:
1345         IF( nn_components /= jp_iam_opa ) THEN
1346            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST)
1347            IF( srcv(jpr_snow  )%laction ) THEN
1348               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean
1349            ENDIF
1350         ENDIF
1351         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:)
1352         ELSE                   ;   qns(:,:) =                              zqns(:,:)
1353         ENDIF
1354
1355         !                                                       ! solar flux over the ocean          (qsr)
1356         IF     ( srcv(jpr_qsroce)%laction ) THEN   ;   zqsr(:,:) = frcv(jpr_qsroce)%z3(:,:,1)
1357         ELSE IF( srcv(jpr_qsrmix)%laction ) then   ;   zqsr(:,:) = frcv(jpr_qsrmix)%z3(:,:,1)
1358         ELSE                                       ;   zqsr(:,:) = 0._wp
1359         ENDIF
1360         IF( ln_dm2dc .AND. ln_cpl )   zqsr(:,:) = sbc_dcy( zqsr )   ! modify qsr to include the diurnal cycle
1361         IF( ln_mixcpl ) THEN   ;   qsr(:,:) = qsr(:,:) * xcplmask(:,:,0) + zqsr(:,:) * zmsk(:,:)
1362         ELSE                   ;   qsr(:,:) =                              zqsr(:,:)
1363         ENDIF
1364         !
1365         ! salt flux over the ocean (received by opa in case of opa <-> sas coupling)
1366         IF( srcv(jpr_sflx )%laction )   sfx(:,:) = frcv(jpr_sflx  )%z3(:,:,1)
1367         ! Ice cover  (received by opa in case of opa <-> sas coupling)
1368         IF( srcv(jpr_fice )%laction )   fr_i(:,:) = frcv(jpr_fice )%z3(:,:,1)
1369         !
1370
1371      ENDIF
1372     
1373      !                                                        ! land ice masses : Greenland
1374      zepsilon = rn_iceshelf_fluxes_tolerance
1375
1376      IF( srcv(jpr_grnm)%laction ) THEN
1377         greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
1378         ! take average over ocean points of input array to avoid cumulative error over time
1379         zgreenland_icesheet_mass_in = SUM( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
1380         IF(lk_mpp) CALL mpp_sum( zgreenland_icesheet_mass_in )
1381         zmask_sum = SUM( tmask(:,:,1) )
1382         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
1383         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum
1384         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
1385         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
1386            zgreenland_icesheet_mass_b = greenland_icesheet_mass
1387           
1388            ! Only update the mass if it has increased
1389            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
1390               greenland_icesheet_mass = zgreenland_icesheet_mass_in
1391            ENDIF
1392           
1393            IF( zgreenland_icesheet_mass_b /= 0.0 ) &
1394           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed 
1395            greenland_icesheet_timelapsed = 0.0_wp       
1396         ENDIF
1397         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
1398         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass
1399         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
1400         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
1401      ENDIF
1402
1403      !                                                        ! land ice masses : Antarctica
1404      IF( srcv(jpr_antm)%laction ) THEN
1405         antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
1406         ! take average over ocean points of input array to avoid cumulative error from rounding errors over time
1407         zantarctica_icesheet_mass_in = SUM( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
1408         IF(lk_mpp) CALL mpp_sum( zantarctica_icesheet_mass_in )
1409         zmask_sum = SUM( tmask(:,:,1) )
1410         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
1411         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum
1412         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
1413         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
1414            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
1415           
1416            ! Only update the mass if it has increased
1417            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
1418               antarctica_icesheet_mass = zantarctica_icesheet_mass_in
1419            END IF
1420           
1421            IF( zantarctica_icesheet_mass_b /= 0.0 ) &
1422          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed 
1423            antarctica_icesheet_timelapsed = 0.0_wp       
1424         ENDIF
1425         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
1426         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass
1427         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
1428         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
1429      ENDIF
1430
1431      !
1432      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2 )
1433      !
1434      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_rcv')
1435      !
1436   END SUBROUTINE sbc_cpl_rcv
1437   
1438
1439   SUBROUTINE sbc_cpl_ice_tau( p_taui, p_tauj )     
1440      !!----------------------------------------------------------------------
1441      !!             ***  ROUTINE sbc_cpl_ice_tau  ***
1442      !!
1443      !! ** Purpose :   provide the stress over sea-ice in coupled mode
1444      !!
1445      !! ** Method  :   transform the received stress from the atmosphere into
1446      !!             an atmosphere-ice stress in the (i,j) ocean referencial
1447      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
1448      !!                'C'-grid : i- (j-) components given at U- (V-) point
1449      !!                'I'-grid : B-grid lower-left corner: both components given at I-point
1450      !!
1451      !!                The received stress are :
1452      !!                 - defined by 3 components (if cartesian coordinate)
1453      !!                        or by 2 components (if spherical)
1454      !!                 - oriented along geographical   coordinate (if eastward-northward)
1455      !!                        or  along the local grid coordinate (if local grid)
1456      !!                 - given at U- and V-point, resp.   if received on 2 grids
1457      !!                        or at a same point (T or I) if received on 1 grid
1458      !!                Therefore and if necessary, they are successively
1459      !!             processed in order to obtain them
1460      !!                 first  as  2 components on the sphere
1461      !!                 second as  2 components oriented along the local grid
1462      !!                 third  as  2 components on the cp_ice_msh point
1463      !!
1464      !!                Except in 'oce and ice' case, only one vector stress field
1465      !!             is received. It has already been processed in sbc_cpl_rcv
1466      !!             so that it is now defined as (i,j) components given at U-
1467      !!             and V-points, respectively. Therefore, only the third
1468      !!             transformation is done and only if the ice-grid is a 'I'-grid.
1469      !!
1470      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
1471      !!----------------------------------------------------------------------
1472      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
1473      REAL(wp), INTENT(out), DIMENSION(:,:) ::   p_tauj   ! at I-point (B-grid) or U & V-point (C-grid)
1474      !!
1475      INTEGER ::   ji, jj                          ! dummy loop indices
1476      INTEGER ::   itx                             ! index of taux over ice
1477      REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty 
1478      !!----------------------------------------------------------------------
1479      !
1480      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_tau')
1481      !
1482      CALL wrk_alloc( jpi,jpj, ztx, zty )
1483
1484      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
1485      ELSE                                ;   itx =  jpr_otx1
1486      ENDIF
1487
1488      ! do something only if we just received the stress from atmosphere
1489      IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
1490
1491         !                                                      ! ======================= !
1492         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
1493            !                                                   ! ======================= !
1494           
1495            IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
1496               !                                                       ! (cartesian to spherical -> 3 to 2 components)
1497               CALL geo2oce(  frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), frcv(jpr_itz1)%z3(:,:,1),   &
1498                  &          srcv(jpr_itx1)%clgrid, ztx, zty )
1499               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 1st grid
1500               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 1st grid
1501               !
1502               IF( srcv(jpr_itx2)%laction ) THEN
1503                  CALL geo2oce( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), frcv(jpr_itz2)%z3(:,:,1),   &
1504                     &          srcv(jpr_itx2)%clgrid, ztx, zty )
1505                  frcv(jpr_itx2)%z3(:,:,1) = ztx(:,:)   ! overwrite 1st comp. on the 2nd grid
1506                  frcv(jpr_ity2)%z3(:,:,1) = zty(:,:)   ! overwrite 2nd comp. on the 2nd grid
1507               ENDIF
1508               !
1509            ENDIF
1510            !
1511            IF( TRIM( sn_rcv_tau%clvor ) == 'eastward-northward' ) THEN   ! 2 components oriented along the local grid
1512               !                                                       ! (geographical to local grid -> rotate the components)
1513               CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->i', ztx )   
1514               IF( srcv(jpr_itx2)%laction ) THEN
1515                  CALL rot_rep( frcv(jpr_itx2)%z3(:,:,1), frcv(jpr_ity2)%z3(:,:,1), srcv(jpr_itx2)%clgrid, 'en->j', zty )   
1516               ELSE
1517                  CALL rot_rep( frcv(jpr_itx1)%z3(:,:,1), frcv(jpr_ity1)%z3(:,:,1), srcv(jpr_itx1)%clgrid, 'en->j', zty ) 
1518               ENDIF
1519               frcv(jpr_itx1)%z3(:,:,1) = ztx(:,:)      ! overwrite 1st component on the 1st grid
1520               frcv(jpr_ity1)%z3(:,:,1) = zty(:,:)      ! overwrite 2nd component on the 1st grid
1521            ENDIF
1522            !                                                   ! ======================= !
1523         ELSE                                                   !     use ocean stress    !
1524            !                                                   ! ======================= !
1525            frcv(jpr_itx1)%z3(:,:,1) = frcv(jpr_otx1)%z3(:,:,1)
1526            frcv(jpr_ity1)%z3(:,:,1) = frcv(jpr_oty1)%z3(:,:,1)
1527            !
1528         ENDIF
1529         !                                                      ! ======================= !
1530         !                                                      !     put on ice grid     !
1531         !                                                      ! ======================= !
1532         !   
1533         !                                                  j+1   j     -----V---F
1534         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
1535         ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
1536         !                                                               |       |
1537         !                                                   j    j-1   -I-------|
1538         !                                               (for I)         |       |
1539         !                                                              i-1  i   i
1540         !                                                               i      i+1 (for I)
1541         SELECT CASE ( cp_ice_msh )
1542            !
1543         CASE( 'I' )                                         ! B-grid ==> I
1544            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1545            CASE( 'U' )
1546               DO jj = 2, jpjm1                                   ! (U,V) ==> I
1547                  DO ji = 2, jpim1   ! NO vector opt.
1548                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji-1,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) )
1549                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1550                  END DO
1551               END DO
1552            CASE( 'F' )
1553               DO jj = 2, jpjm1                                   ! F ==> I
1554                  DO ji = 2, jpim1   ! NO vector opt.
1555                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji-1,jj-1,1)
1556                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji-1,jj-1,1)
1557                  END DO
1558               END DO
1559            CASE( 'T' )
1560               DO jj = 2, jpjm1                                   ! T ==> I
1561                  DO ji = 2, jpim1   ! NO vector opt.
1562                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji-1,jj  ,1)   &
1563                        &                   + frcv(jpr_itx1)%z3(ji,jj-1,1) + frcv(jpr_itx1)%z3(ji-1,jj-1,1) ) 
1564                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1)   &
1565                        &                   + frcv(jpr_oty1)%z3(ji,jj-1,1) + frcv(jpr_ity1)%z3(ji-1,jj-1,1) )
1566                  END DO
1567               END DO
1568            CASE( 'I' )
1569               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! I ==> I
1570               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1571            END SELECT
1572            IF( srcv(jpr_itx1)%clgrid /= 'I' ) THEN
1573               CALL lbc_lnk( p_taui, 'I',  -1. )   ;   CALL lbc_lnk( p_tauj, 'I',  -1. )
1574            ENDIF
1575            !
1576         CASE( 'F' )                                         ! B-grid ==> F
1577            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1578            CASE( 'U' )
1579               DO jj = 2, jpjm1                                   ! (U,V) ==> F
1580                  DO ji = fs_2, fs_jpim1   ! vector opt.
1581                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj+1,1) )
1582                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1) )
1583                  END DO
1584               END DO
1585            CASE( 'I' )
1586               DO jj = 2, jpjm1                                   ! I ==> F
1587                  DO ji = 2, jpim1   ! NO vector opt.
1588                     p_taui(ji,jj) = frcv(jpr_itx1)%z3(ji+1,jj+1,1)
1589                     p_tauj(ji,jj) = frcv(jpr_ity1)%z3(ji+1,jj+1,1)
1590                  END DO
1591               END DO
1592            CASE( 'T' )
1593               DO jj = 2, jpjm1                                   ! T ==> F
1594                  DO ji = 2, jpim1   ! NO vector opt.
1595                     p_taui(ji,jj) = 0.25 * ( frcv(jpr_itx1)%z3(ji,jj  ,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1)   &
1596                        &                   + frcv(jpr_itx1)%z3(ji,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj+1,1) ) 
1597                     p_tauj(ji,jj) = 0.25 * ( frcv(jpr_ity1)%z3(ji,jj  ,1) + frcv(jpr_ity1)%z3(ji+1,jj  ,1)   &
1598                        &                   + frcv(jpr_ity1)%z3(ji,jj+1,1) + frcv(jpr_ity1)%z3(ji+1,jj+1,1) )
1599                  END DO
1600               END DO
1601            CASE( 'F' )
1602               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! F ==> F
1603               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1604            END SELECT
1605            IF( srcv(jpr_itx1)%clgrid /= 'F' ) THEN
1606               CALL lbc_lnk( p_taui, 'F',  -1. )   ;   CALL lbc_lnk( p_tauj, 'F',  -1. )
1607            ENDIF
1608            !
1609         CASE( 'C' )                                         ! C-grid ==> U,V
1610            SELECT CASE ( srcv(jpr_itx1)%clgrid )
1611            CASE( 'U' )
1612               p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
1613               p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
1614            CASE( 'F' )
1615               DO jj = 2, jpjm1                                   ! F ==> (U,V)
1616                  DO ji = fs_2, fs_jpim1   ! vector opt.
1617                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
1618                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(jj,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
1619                  END DO
1620               END DO
1621            CASE( 'T' )
1622               DO jj = 2, jpjm1                                   ! T ==> (U,V)
1623                  DO ji = fs_2, fs_jpim1   ! vector opt.
1624                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
1625                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
1626                  END DO
1627               END DO
1628            CASE( 'I' )
1629               DO jj = 2, jpjm1                                   ! I ==> (U,V)
1630                  DO ji = 2, jpim1   ! NO vector opt.
1631                     p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
1632                     p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
1633                  END DO
1634               END DO
1635            END SELECT
1636            IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN
1637               CALL lbc_lnk( p_taui, 'U',  -1. )   ;   CALL lbc_lnk( p_tauj, 'V',  -1. )
1638            ENDIF
1639         END SELECT
1640
1641      ENDIF
1642      !   
1643      CALL wrk_dealloc( jpi,jpj, ztx, zty )
1644      !
1645      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_tau')
1646      !
1647   END SUBROUTINE sbc_cpl_ice_tau
1648   
1649
1650   SUBROUTINE sbc_cpl_ice_flx( p_frld, palbi, psst, pist )
1651      !!----------------------------------------------------------------------
1652      !!             ***  ROUTINE sbc_cpl_ice_flx  ***
1653      !!
1654      !! ** Purpose :   provide the heat and freshwater fluxes of the
1655      !!              ocean-ice system.
1656      !!
1657      !! ** Method  :   transform the fields received from the atmosphere into
1658      !!             surface heat and fresh water boundary condition for the
1659      !!             ice-ocean system. The following fields are provided:
1660      !!              * total non solar, solar and freshwater fluxes (qns_tot,
1661      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
1662      !!             NB: emp_tot include runoffs and calving.
1663      !!              * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
1664      !!             emp_ice = sublimation - solid precipitation as liquid
1665      !!             precipitation are re-routed directly to the ocean and
1666      !!             runoffs and calving directly enter the ocean.
1667      !!              * solid precipitation (sprecip), used to add to qns_tot
1668      !!             the heat lost associated to melting solid precipitation
1669      !!             over the ocean fraction.
1670      !!       ===>> CAUTION here this changes the net heat flux received from
1671      !!             the atmosphere
1672      !!
1673      !!                  - the fluxes have been separated from the stress as
1674      !!                 (a) they are updated at each ice time step compare to
1675      !!                 an update at each coupled time step for the stress, and
1676      !!                 (b) the conservative computation of the fluxes over the
1677      !!                 sea-ice area requires the knowledge of the ice fraction
1678      !!                 after the ice advection and before the ice thermodynamics,
1679      !!                 so that the stress is updated before the ice dynamics
1680      !!                 while the fluxes are updated after it.
1681      !!
1682      !! ** Action  :   update at each nf_ice time step:
1683      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
1684      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
1685      !!                   emp_tot            total evaporation - precipitation(liquid and solid) (-runoff)(-calving)
1686      !!                   emp_ice            ice sublimation - solid precipitation over the ice
1687      !!                   dqns_ice           d(non-solar heat flux)/d(Temperature) over the ice
1688      !!                   sprecip             solid precipitation over the ocean 
1689      !!----------------------------------------------------------------------
1690      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   p_frld     ! lead fraction                [0 to 1]
1691      ! optional arguments, used only in 'mixed oce-ice' case
1692      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo
1693      REAL(wp), INTENT(in   ), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
1694      REAL(wp), INTENT(in   ), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin]
1695      !
1696      INTEGER ::   jl         ! dummy loop index
1697      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, ztmp, zicefr, zmsk
1698      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot
1699      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice
1700      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zevap, zsnw, zqns_oce, zqsr_oce, zqprec_ice, zqemp_oce ! for LIM3
1701      !!----------------------------------------------------------------------
1702      !
1703      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx')
1704      !
1705      CALL wrk_alloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
1706      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
1707
1708      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
1709      zicefr(:,:) = 1.- p_frld(:,:)
1710      zcptn(:,:) = rcp * sst_m(:,:)
1711      !
1712      !                                                      ! ========================= !
1713      !                                                      !    freshwater budget      !   (emp)
1714      !                                                      ! ========================= !
1715      !
1716      !                                                           ! total Precipitation - total Evaporation (emp_tot)
1717      !                                                           ! solid precipitation - sublimation       (emp_ice)
1718      !                                                           ! solid Precipitation                     (sprecip)
1719      !                                                           ! liquid + solid Precipitation            (tprecip)
1720      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
1721      CASE( 'conservative'  )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
1722         zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
1723         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
1724         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)         
1725#if defined key_cice
1726         IF ( TRIM(sn_rcv_emp%clcat) == 'yes' ) THEN
1727            ! zemp_ice is the sum of frcv(jpr_ievp)%z3(:,:,1) over all layers - snow
1728            zemp_ice(:,:) = - frcv(jpr_snow)%z3(:,:,1)
1729            DO jl=1,jpl
1730               zemp_ice(:,:   ) = zemp_ice(:,:) + frcv(jpr_ievp)%z3(:,:,jl)
1731            ENDDO
1732            ! latent heat coupled for each category in CICE
1733            qla_ice(:,:,1:jpl) = - frcv(jpr_ievp)%z3(:,:,1:jpl) * lsub
1734         ELSE
1735            ! If CICE has multicategories it still expects coupling fields for
1736            ! each even if we treat as a single field
1737            ! The latent heat flux is split between the ice categories according
1738            ! to the fraction of the ice in each category
1739            zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1740            WHERE ( zicefr(:,:) /= 0._wp ) 
1741               ztmp(:,:) = 1./zicefr(:,:)
1742            ELSEWHERE 
1743               ztmp(:,:) = 0.e0
1744            END WHERE 
1745            DO jl=1,jpl
1746               qla_ice(:,:,jl) = - a_i(:,:,jl) * ztmp(:,:) * frcv(jpr_ievp)%z3(:,:,1) * lsub 
1747            END DO
1748            WHERE ( zicefr(:,:) == 0._wp )  qla_ice(:,:,1) = -frcv(jpr_ievp)%z3(:,:,1) * lsub 
1749         ENDIF
1750#else         
1751         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1)
1752#endif                 
1753            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation
1754         IF( iom_use('hflx_rain_cea') )   &
1755            CALL iom_put( 'hflx_rain_cea', frcv(jpr_rain)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from liq. precip.
1756         IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') )   &
1757            ztmp(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:)
1758         IF( iom_use('evap_ao_cea'  ) )   &
1759            CALL iom_put( 'evap_ao_cea'  , ztmp                   )   ! ice-free oce evap (cell average)
1760         IF( iom_use('hflx_evap_cea') )   &
1761            CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * zcptn(:,:) )   ! heat flux from from evap (cell average)
1762      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
1763         zemp_tot(:,:) = p_frld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + zicefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
1764         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1)
1765         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
1766         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
1767      END SELECT
1768
1769      IF( iom_use('subl_ai_cea') )   &
1770         CALL iom_put( 'subl_ai_cea', frcv(jpr_ievp)%z3(:,:,1) * zicefr(:,:) )   ! Sublimation over sea-ice         (cell average)
1771      !   
1772      !                                                           ! runoffs and calving (put in emp_tot)
1773      IF( srcv(jpr_rnf)%laction )   rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
1774      IF( srcv(jpr_cal)%laction ) THEN
1775         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
1776         CALL iom_put( 'calving_cea', frcv(jpr_cal)%z3(:,:,1) )
1777      ENDIF
1778
1779      IF( ln_mixcpl ) THEN
1780         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
1781         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
1782         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
1783         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
1784      ELSE
1785         emp_tot(:,:) =                                  zemp_tot(:,:)
1786         emp_ice(:,:) =                                  zemp_ice(:,:)
1787         sprecip(:,:) =                                  zsprecip(:,:)
1788         tprecip(:,:) =                                  ztprecip(:,:)
1789      ENDIF
1790
1791         CALL iom_put( 'snowpre'    , sprecip                                )   ! Snow
1792      IF( iom_use('snow_ao_cea') )   &
1793         CALL iom_put( 'snow_ao_cea', sprecip(:,:) * p_frld(:,:)             )   ! Snow        over ice-free ocean  (cell average)
1794      IF( iom_use('snow_ai_cea') )   &
1795         CALL iom_put( 'snow_ai_cea', sprecip(:,:) * zicefr(:,:)             )   ! Snow        over sea-ice         (cell average)
1796
1797      !                                                      ! ========================= !
1798      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
1799      !                                                      ! ========================= !
1800      CASE( 'oce only' )                                     ! the required field is directly provided
1801         zqns_tot(:,:  ) = frcv(jpr_qnsoce)%z3(:,:,1)
1802      CASE( 'conservative' )                                      ! the required fields are directly provided
1803         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1804         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1805            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
1806         ELSE
1807            ! Set all category values equal for the moment
1808            DO jl=1,jpl
1809               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1810            ENDDO
1811         ENDIF
1812      CASE( 'oce and ice' )       ! the total flux is computed from ocean and ice fluxes
1813         zqns_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
1814         IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
1815            DO jl=1,jpl
1816               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)   
1817               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
1818            ENDDO
1819         ELSE
1820            qns_tot(:,:   ) = qns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1821            DO jl=1,jpl
1822               zqns_tot(:,:   ) = zqns_tot(:,:) + zicefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
1823               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
1824            ENDDO
1825         ENDIF
1826      CASE( 'mixed oce-ice' )     ! the ice flux is cumputed from the total flux, the SST and ice informations
1827! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1828         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
1829         zqns_ice(:,:,1) = frcv(jpr_qnsmix)%z3(:,:,1)    &
1830            &            + frcv(jpr_dqnsdt)%z3(:,:,1) * ( pist(:,:,1) - ( (rt0 + psst(:,:  ) ) * p_frld(:,:)   &
1831            &                                                   +          pist(:,:,1)   * zicefr(:,:) ) )
1832      END SELECT
1833!!gm
1834!!    currently it is taken into account in leads budget but not in the zqns_tot, and thus not in
1835!!    the flux that enter the ocean....
1836!!    moreover 1 - it is not diagnose anywhere....
1837!!             2 - it is unclear for me whether this heat lost is taken into account in the atmosphere or not...
1838!!
1839!! similar job should be done for snow and precipitation temperature
1840      !                                     
1841      IF( srcv(jpr_cal)%laction ) THEN                            ! Iceberg melting
1842         ztmp(:,:) = frcv(jpr_cal)%z3(:,:,1) * lfus               ! add the latent heat of iceberg melting
1843         zqns_tot(:,:) = zqns_tot(:,:) - ztmp(:,:)
1844         IF( iom_use('hflx_cal_cea') )   &
1845            CALL iom_put( 'hflx_cal_cea', ztmp + frcv(jpr_cal)%z3(:,:,1) * zcptn(:,:) )   ! heat flux from calving
1846      ENDIF
1847
1848      ztmp(:,:) = p_frld(:,:) * zsprecip(:,:) * lfus
1849      IF( iom_use('hflx_snow_cea') )    CALL iom_put( 'hflx_snow_cea', ztmp + sprecip(:,:) * zcptn(:,:) )   ! heat flux from snow (cell average)
1850
1851#if defined key_lim3
1852      CALL wrk_alloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
1853
1854      ! --- evaporation --- !
1855      ! clem: evap_ice is set to 0 for LIM3 since we still do not know what to do with sublimation
1856      ! the problem is: the atm. imposes both mass evaporation and heat removed from the snow/ice
1857      !                 but it is incoherent WITH the ice model 
1858      DO jl=1,jpl
1859         evap_ice(:,:,jl) = 0._wp  ! should be: frcv(jpr_ievp)%z3(:,:,1)
1860      ENDDO
1861      zevap(:,:) = zemp_tot(:,:) + ztprecip(:,:) ! evaporation over ocean
1862
1863      ! --- evaporation minus precipitation --- !
1864      emp_oce(:,:) = emp_tot(:,:) - emp_ice(:,:)
1865
1866      ! --- non solar flux over ocean --- !
1867      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1868      zqns_oce = 0._wp
1869      WHERE( p_frld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / p_frld(:,:)
1870
1871      ! --- heat flux associated with emp --- !
1872      zsnw(:,:) = 0._wp
1873      CALL lim_thd_snwblow( p_frld, zsnw )  ! snow distribution over ice after wind blowing
1874      zqemp_oce(:,:) = -      zevap(:,:)                   * p_frld(:,:)      *   zcptn(:,:)   &      ! evap
1875         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptn(:,:)   &      ! liquid precip
1876         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptn(:,:) - lfus ) ! solid precip over ocean
1877      qemp_ice(:,:)  = -   frcv(jpr_ievp)%z3(:,:,1)        * zicefr(:,:)      *   zcptn(:,:)   &      ! ice evap
1878         &             +   zsprecip(:,:)                   * zsnw             * ( zcptn(:,:) - lfus ) ! solid precip over ice
1879
1880      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1881      zqprec_ice(:,:) = rhosn * ( zcptn(:,:) - lfus )
1882
1883      ! --- total non solar flux --- !
1884      zqns_tot(:,:) = zqns_tot(:,:) + qemp_ice(:,:) + zqemp_oce(:,:)
1885
1886      ! --- in case both coupled/forced are active, we must mix values --- !
1887      IF( ln_mixcpl ) THEN
1888         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
1889         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
1890         DO jl=1,jpl
1891            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1892         ENDDO
1893         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
1894         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
1895!!clem         evap_ice(:,:) = evap_ice(:,:) * xcplmask(:,:,0)
1896      ELSE
1897         qns_tot  (:,:  ) = zqns_tot  (:,:  )
1898         qns_oce  (:,:  ) = zqns_oce  (:,:  )
1899         qns_ice  (:,:,:) = zqns_ice  (:,:,:)
1900         qprec_ice(:,:)   = zqprec_ice(:,:)
1901         qemp_oce (:,:)   = zqemp_oce (:,:)
1902      ENDIF
1903
1904      CALL wrk_dealloc( jpi,jpj, zevap, zsnw, zqns_oce, zqprec_ice, zqemp_oce ) 
1905#else
1906
1907      ! clem: this formulation is certainly wrong... but better than it was...
1908      zqns_tot(:,:) = zqns_tot(:,:)                       &            ! zqns_tot update over free ocean with:
1909         &          - ztmp(:,:)                           &            ! remove the latent heat flux of solid precip. melting
1910         &          - (  zemp_tot(:,:)                    &            ! remove the heat content of mass flux (assumed to be at SST)
1911         &             - zemp_ice(:,:) * zicefr(:,:)  ) * zcptn(:,:) 
1912
1913     IF( ln_mixcpl ) THEN
1914         qns_tot(:,:) = qns(:,:) * p_frld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1915         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
1916         DO jl=1,jpl
1917            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
1918         ENDDO
1919      ELSE
1920         qns_tot(:,:  ) = zqns_tot(:,:  )
1921         qns_ice(:,:,:) = zqns_ice(:,:,:)
1922      ENDIF
1923
1924#endif
1925
1926      !                                                      ! ========================= !
1927      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
1928      !                                                      ! ========================= !
1929      CASE( 'oce only' )
1930         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
1931      CASE( 'conservative' )
1932         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1933         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1934            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
1935         ELSE
1936            ! Set all category values equal for the moment
1937            DO jl=1,jpl
1938               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1939            ENDDO
1940         ENDIF
1941         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1942         zqsr_ice(:,:,1) = frcv(jpr_qsrice)%z3(:,:,1)
1943      CASE( 'oce and ice' )
1944         zqsr_tot(:,:  ) =  p_frld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
1945         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
1946            DO jl=1,jpl
1947               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)   
1948               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
1949            ENDDO
1950         ELSE
1951            qsr_tot(:,:   ) = qsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1952            DO jl=1,jpl
1953               zqsr_tot(:,:   ) = zqsr_tot(:,:) + zicefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
1954               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
1955            ENDDO
1956         ENDIF
1957      CASE( 'mixed oce-ice' )
1958         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
1959! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
1960!       Create solar heat flux over ice using incoming solar heat flux and albedos
1961!       ( see OASIS3 user guide, 5th edition, p39 )
1962         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   &
1963            &            / (  1.- ( albedo_oce_mix(:,:  ) * p_frld(:,:)       &
1964            &                     + palbi         (:,:,1) * zicefr(:,:) ) )
1965      END SELECT
1966      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
1967         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
1968         DO jl=1,jpl
1969            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
1970         ENDDO
1971      ENDIF
1972
1973#if defined key_lim3
1974      CALL wrk_alloc( jpi,jpj, zqsr_oce ) 
1975      ! --- solar flux over ocean --- !
1976      !         note: p_frld cannot be = 0 since we limit the ice concentration to amax
1977      zqsr_oce = 0._wp
1978      WHERE( p_frld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / p_frld(:,:)
1979
1980      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
1981      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF
1982
1983      CALL wrk_dealloc( jpi,jpj, zqsr_oce ) 
1984#endif
1985
1986      IF( ln_mixcpl ) THEN
1987         qsr_tot(:,:) = qsr(:,:) * p_frld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
1988         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
1989         DO jl=1,jpl
1990            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
1991         ENDDO
1992      ELSE
1993         qsr_tot(:,:  ) = zqsr_tot(:,:  )
1994         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
1995      ENDIF
1996
1997      !                                                      ! ========================= !
1998      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
1999      !                                                      ! ========================= !
2000      CASE ('coupled')
2001         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
2002            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
2003         ELSE
2004            ! Set all category values equal for the moment
2005            DO jl=1,jpl
2006               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
2007            ENDDO
2008         ENDIF
2009      END SELECT
2010     
2011      IF( ln_mixcpl ) THEN
2012         DO jl=1,jpl
2013            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
2014         ENDDO
2015      ELSE
2016         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
2017      ENDIF
2018     
2019      !                                                      ! ========================= !
2020      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !    topmelt and botmelt    !
2021      !                                                      ! ========================= !
2022      CASE ('coupled')
2023         topmelt(:,:,:)=frcv(jpr_topm)%z3(:,:,:)
2024         botmelt(:,:,:)=frcv(jpr_botm)%z3(:,:,:)
2025      END SELECT
2026
2027      ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 )
2028      ! Used for LIM2 and LIM3
2029      ! Coupled case: since cloud cover is not received from atmosphere
2030      !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
2031      fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )
2032      fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )
2033
2034      CALL wrk_dealloc( jpi,jpj,     zcptn, ztmp, zicefr, zmsk, zemp_tot, zemp_ice, zsprecip, ztprecip, zqns_tot, zqsr_tot )
2035      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice )
2036      !
2037      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_ice_flx')
2038      !
2039   END SUBROUTINE sbc_cpl_ice_flx
2040   
2041   
2042   SUBROUTINE sbc_cpl_snd( kt )
2043      !!----------------------------------------------------------------------
2044      !!             ***  ROUTINE sbc_cpl_snd  ***
2045      !!
2046      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
2047      !!
2048      !! ** Method  :   send to the atmosphere through a call to cpl_snd
2049      !!              all the needed fields (as defined in sbc_cpl_init)
2050      !!----------------------------------------------------------------------
2051      INTEGER, INTENT(in) ::   kt
2052      !
2053      INTEGER ::   ji, jj, jl   ! dummy loop indices
2054      INTEGER ::   ikchoix
2055      INTEGER ::   isec, info   ! local integer
2056      REAL(wp) ::   zumax, zvmax
2057      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
2058      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
2059
2060#if defined key_medusa 
2061! RSRH Temporarily definitition of array used in medusa coupling until we actually have the medusa code included in the build
2062      REAL(wp)  ::  trc2d_TEST(jpi,jpj) 
2063#endif 
2064      !!----------------------------------------------------------------------
2065      !
2066      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_snd')
2067      !
2068      CALL wrk_alloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2069      CALL wrk_alloc( jpi,jpj,jpl, ztmp3, ztmp4 )
2070
2071      isec = ( kt - nit000 ) * NINT(rdttra(1))        ! date of exchanges
2072
2073      zfr_l(:,:) = 1.- fr_i(:,:)
2074      !                                                      ! ------------------------- !
2075      !                                                      !    Surface temperature    !   in Kelvin
2076      !                                                      ! ------------------------- !
2077      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN
2078         
2079         IF ( nn_components == jp_iam_opa ) THEN
2080            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of ln_useCT on the received part
2081         ELSE
2082            ! we must send the surface potential temperature
2083            IF( ln_useCT )  THEN    ;   ztmp1(:,:) = eos_pt_from_ct( tsn(:,:,1,jp_tem), tsn(:,:,1,jp_sal) )
2084            ELSE                    ;   ztmp1(:,:) = tsn(:,:,1,jp_tem)
2085            ENDIF
2086            !
2087            SELECT CASE( sn_snd_temp%cldes)
2088            CASE( 'oce only'             )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2089            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
2090               SELECT CASE( sn_snd_temp%clcat )
2091               CASE( 'yes' )   
2092                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
2093               CASE( 'no' )
2094                  WHERE( SUM( a_i, dim=3 ) /= 0. )
2095                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
2096                  ELSEWHERE
2097                     ztmp3(:,:,1) = rt0 ! TODO: Is freezing point a good default? (Maybe SST is better?)
2098                  END WHERE
2099               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2100               END SELECT
2101            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)   
2102               SELECT CASE( sn_snd_temp%clcat )
2103               CASE( 'yes' )   
2104                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2105               CASE( 'no' )
2106                  ztmp3(:,:,:) = 0.0
2107                  DO jl=1,jpl
2108                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
2109                  ENDDO
2110               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2111               END SELECT
2112            CASE( 'oce and weighted ice' )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0 
2113               SELECT CASE( sn_snd_temp%clcat )
2114               CASE( 'yes' )   
2115           ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2116               CASE( 'no' )
2117           ztmp3(:,:,:) = 0.0
2118           DO jl=1,jpl
2119                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
2120           ENDDO
2121               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
2122               END SELECT
2123            CASE( 'mixed oce-ice'        )   
2124               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:) 
2125               DO jl=1,jpl
2126                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
2127               ENDDO
2128            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
2129            END SELECT
2130         ENDIF
2131         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2132         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
2133         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2134      ENDIF
2135      !                                                      ! ------------------------- !
2136      !                                                      !           Albedo          !
2137      !                                                      ! ------------------------- !
2138      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
2139         SELECT CASE( sn_snd_alb%cldes )
2140         CASE( 'ice'          )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
2141         CASE( 'weighted ice' )   ; ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2142         CASE default             ; CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
2143         END SELECT
2144         CALL cpl_snd( jps_albice, isec, ztmp3, info )
2145      ENDIF
2146      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
2147         ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:)
2148         DO jl=1,jpl
2149            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
2150         ENDDO
2151         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2152      ENDIF
2153      !                                                      ! ------------------------- !
2154      !                                                      !  Ice fraction & Thickness !
2155      !                                                      ! ------------------------- !
2156      ! Send ice fraction field to atmosphere
2157      IF( ssnd(jps_fice)%laction ) THEN
2158         SELECT CASE( sn_snd_thick%clcat )
2159         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
2160         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
2161         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2162         END SELECT
2163         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info )
2164      ENDIF
2165     
2166      ! Send ice fraction field (first order interpolation), for weighting UM fluxes to be passed to NEMO
2167      IF (ssnd(jps_fice1)%laction) THEN
2168         SELECT CASE (sn_snd_thick1%clcat)
2169         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
2170         CASE( 'no' )    ;   ztmp3(:,:,1) = fr_i(:,:)
2171         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' )
2172    END SELECT
2173         CALL cpl_snd (jps_fice1, isec, ztmp3, info)
2174      ENDIF
2175     
2176      ! Send ice fraction field to OPA (sent by SAS in SAS-OPA coupling)
2177      IF( ssnd(jps_fice2)%laction ) THEN
2178         ztmp3(:,:,1) = fr_i(:,:)
2179         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
2180      ENDIF
2181
2182      ! Send ice and snow thickness field
2183      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
2184         SELECT CASE( sn_snd_thick%cldes)
2185         CASE( 'none'                  )       ! nothing to do
2186         CASE( 'weighted ice and snow' )   
2187            SELECT CASE( sn_snd_thick%clcat )
2188            CASE( 'yes' )   
2189               ztmp3(:,:,1:jpl) =  ht_i(:,:,1:jpl) * a_i(:,:,1:jpl)
2190               ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
2191            CASE( 'no' )
2192               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
2193               DO jl=1,jpl
2194                  ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
2195                  ztmp4(:,:,1) = ztmp4(:,:,1) + ht_s(:,:,jl) * a_i(:,:,jl)
2196               ENDDO
2197            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2198            END SELECT
2199         CASE( 'ice and snow'         )   
2200            SELECT CASE( sn_snd_thick%clcat )
2201            CASE( 'yes' )
2202               ztmp3(:,:,1:jpl) = ht_i(:,:,1:jpl)
2203               ztmp4(:,:,1:jpl) = ht_s(:,:,1:jpl)
2204            CASE( 'no' )
2205               WHERE( SUM( a_i, dim=3 ) /= 0. )
2206                  ztmp3(:,:,1) = SUM( ht_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
2207                  ztmp4(:,:,1) = SUM( ht_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
2208               ELSEWHERE
2209                 ztmp3(:,:,1) = 0.
2210                 ztmp4(:,:,1) = 0.
2211               END WHERE
2212            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
2213            END SELECT
2214         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
2215         END SELECT
2216         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
2217         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
2218      ENDIF
2219      !
2220      ! Send meltpond fields
2221      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN
2222         SELECT CASE( sn_snd_mpnd%cldes) 
2223         CASE( 'weighted ice' ) 
2224            SELECT CASE( sn_snd_mpnd%clcat ) 
2225            CASE( 'yes' ) 
2226               ztmp3(:,:,1:jpl) =  a_p(:,:,1:jpl) * a_i(:,:,1:jpl) 
2227               ztmp4(:,:,1:jpl) =  ht_p(:,:,1:jpl) * a_i(:,:,1:jpl) 
2228            CASE( 'no' ) 
2229               ztmp3(:,:,:) = 0.0 
2230               ztmp4(:,:,:) = 0.0 
2231               DO jl=1,jpl 
2232                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_p(:,:,jpl) * a_i(:,:,jpl) 
2233                 ztmp4(:,:,1) = ztmp4(:,:,1) + ht_p(:,:,jpl) * a_i(:,:,jpl) 
2234               ENDDO 
2235            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' ) 
2236            END SELECT
2237         CASE( 'ice only' )   
2238            ztmp3(:,:,1:jpl) = a_p(:,:,1:jpl) 
2239            ztmp4(:,:,1:jpl) = ht_p(:,:,1:jpl) 
2240         END SELECT
2241         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )   
2242         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )   
2243         !
2244         ! Send ice effective conductivity
2245         SELECT CASE( sn_snd_cond%cldes)
2246         CASE( 'weighted ice' )   
2247            SELECT CASE( sn_snd_cond%clcat )
2248            CASE( 'yes' )   
2249               ztmp3(:,:,1:jpl) =  kn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
2250            CASE( 'no' )
2251               ztmp3(:,:,:) = 0.0
2252               DO jl=1,jpl
2253                 ztmp3(:,:,1) = ztmp3(:,:,1) + kn_ice(:,:,jl) * a_i(:,:,jl)
2254               ENDDO
2255            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' )
2256            END SELECT
2257         CASE( 'ice only' )   
2258           ztmp3(:,:,1:jpl) = kn_ice(:,:,1:jpl)
2259         END SELECT
2260         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info )
2261      ENDIF
2262      !
2263      !
2264#if defined key_cpl_carbon_cycle
2265      !                                                      ! ------------------------- !
2266      !                                                      !  CO2 flux from PISCES     !
2267      !                                                      ! ------------------------- !
2268      IF( ssnd(jps_co2)%laction )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
2269      !
2270#endif
2271
2272
2273!! JPALM : 03-feb-2015 coupling MEDUSA-UKCA
2274!! add soing exactly the same that CO2 fluxes
2275!! of PISCES
2276!! add CO2-MEDUSA
2277!! add DMS-MEDUSA
2278!! May add also a coupling MED-UKCA key
2279
2280! RSRH set up a test field purely for use during development:
2281      trc2d_TEST(:,:) = 0.0 
2282
2283      ! RSRH. We don't want to use magic numbers in the final code (i.e. 98 and 221).
2284      ! These need moving to a parameter statement (as part of MEDUSA code) or specifying in a namelist
2285      ! so the following code MUST NOT be viewed as anything more than temporary.
2286      IF ( ssnd(jps_bio_co2)%laction ) THEN
2287!         CALL cpl_prism_snd( jps_bio_co2, isec, trc2d(:,:,98:98), info )
2288! RSRH Send our test field in place of the real bio_co2
2289         CALL cpl_prism_snd( jps_bio_co2, isec, trc2d_TEST(:,:), info )
2290      END IF
2291
2292      IF ( ssnd(jps_bio_dms)%laction ) THEN
2293          ! We need to multiply DMS by a conversion factor to get values in the standard units expected in
2294          ! the coupling space.
2295!          ztmp1(:,: ) = trc2d(:,:,221) * dms_unit_conv
2296!RSRH Use our test field in place of the real bio_dms
2297         ztmp1(:,: ) = trc2d_TEST(:,:) * dms_unit_conv
2298         CALL cpl_prism_snd( jps_bio_dms, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2299      ENDIF
2300
2301      !                                                      ! ------------------------- !
2302      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
2303         !                                                   ! ------------------------- !
2304         !   
2305         !                                                  j+1   j     -----V---F
2306         ! surface velocity always sent from T point                     !       |
2307         ! [except for HadGEM3]                                   j      |   T   U
2308         !                                                               |       |
2309         !                                                   j    j-1   -I-------|
2310         !                                               (for I)         |       |
2311         !                                                              i-1  i   i
2312         !                                                               i      i+1 (for I)
2313         IF( nn_components == jp_iam_opa ) THEN
2314            zotx1(:,:) = un(:,:,1) 
2315            zoty1(:,:) = vn(:,:,1) 
2316         ELSE       
2317            SELECT CASE( TRIM( sn_snd_crt%cldes ) )
2318            CASE( 'oce only'             )      ! C-grid ==> T
2319               IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
2320                  DO jj = 2, jpjm1
2321                     DO ji = fs_2, fs_jpim1   ! vector opt.
2322                        zotx1(ji,jj) = 0.5 * ( un(ji,jj,1) + un(ji-1,jj  ,1) )
2323                        zoty1(ji,jj) = 0.5 * ( vn(ji,jj,1) + vn(ji  ,jj-1,1) ) 
2324                     END DO
2325                  END DO
2326               ELSE
2327! Temporarily Changed for UKV
2328                  DO jj = 2, jpjm1
2329                     DO ji = 2, jpim1
2330                        zotx1(ji,jj) = un(ji,jj,1)
2331                        zoty1(ji,jj) = vn(ji,jj,1)
2332                     END DO
2333                  END DO
2334               ENDIF
2335            CASE( 'weighted oce and ice' )   
2336               SELECT CASE ( cp_ice_msh )
2337               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2338                  DO jj = 2, jpjm1
2339                     DO ji = fs_2, fs_jpim1   ! vector opt.
2340                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2341                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)
2342                        zitx1(ji,jj) = 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2343                        zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2344                     END DO
2345                  END DO
2346               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2347                  DO jj = 2, jpjm1
2348                     DO ji = 2, jpim1   ! NO vector opt.
2349                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2350                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2351                        zitx1(ji,jj) = 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2352                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2353                        zity1(ji,jj) = 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2354                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2355                     END DO
2356                  END DO
2357               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2358                  DO jj = 2, jpjm1
2359                     DO ji = 2, jpim1   ! NO vector opt.
2360                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj) 
2361                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj) 
2362                        zitx1(ji,jj) = 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2363                           &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2364                        zity1(ji,jj) = 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2365                           &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2366                     END DO
2367                  END DO
2368               END SELECT
2369               CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
2370            CASE( 'mixed oce-ice'        )
2371               SELECT CASE ( cp_ice_msh )
2372               CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
2373                  DO jj = 2, jpjm1
2374                     DO ji = fs_2, fs_jpim1   ! vector opt.
2375                        zotx1(ji,jj) = 0.5 * ( un   (ji,jj,1) + un   (ji-1,jj  ,1) ) * zfr_l(ji,jj)   &
2376                           &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)
2377                        zoty1(ji,jj) = 0.5 * ( vn   (ji,jj,1) + vn   (ji  ,jj-1,1) ) * zfr_l(ji,jj)   &
2378                           &         + 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)
2379                     END DO
2380                  END DO
2381               CASE( 'I' )                      ! Ocean on C grid, Ice on I-point (B-grid) ==> T
2382                  DO jj = 2, jpjm1
2383                     DO ji = 2, jpim1   ! NO vector opt.
2384                        zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj  ,1) ) * zfr_l(ji,jj)   &   
2385                           &         + 0.25 * ( u_ice(ji+1,jj+1) + u_ice(ji,jj+1)                     &
2386                           &                  + u_ice(ji+1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2387                        zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji  ,jj-1,1) ) * zfr_l(ji,jj)   & 
2388                           &         + 0.25 * ( v_ice(ji+1,jj+1) + v_ice(ji,jj+1)                     &
2389                           &                  + v_ice(ji+1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2390                     END DO
2391                  END DO
2392               CASE( 'F' )                      ! Ocean on C grid, Ice on F-point (B-grid) ==> T
2393                  IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
2394                     DO jj = 2, jpjm1
2395                        DO ji = 2, jpim1   ! NO vector opt.
2396                           zotx1(ji,jj) = 0.5  * ( un(ji,jj,1)      + un(ji-1,jj,1) ) * zfr_l(ji,jj)   &   
2397                                &         + 0.25 * ( u_ice(ji-1,jj-1) + u_ice(ji,jj-1)                     &
2398                                &                  + u_ice(ji-1,jj  ) + u_ice(ji,jj  )  ) *  fr_i(ji,jj)
2399                           zoty1(ji,jj) = 0.5  * ( vn(ji,jj,1)      + vn(ji,jj-1,1) ) * zfr_l(ji,jj)   &
2400                                &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
2401                                &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
2402                        END DO
2403                     END DO
2404#if defined key_cice
2405                  ELSE
2406! Temporarily Changed for HadGEM3
2407                     DO jj = 2, jpjm1
2408                        DO ji = 2, jpim1   ! NO vector opt.
2409                           zotx1(ji,jj) = (1.0-fr_iu(ji,jj)) * un(ji,jj,1)             &
2410                                &              + fr_iu(ji,jj) * 0.5 * ( u_ice(ji,jj-1) + u_ice(ji,jj) ) 
2411                           zoty1(ji,jj) = (1.0-fr_iv(ji,jj)) * vn(ji,jj,1)             &
2412                                &              + fr_iv(ji,jj) * 0.5 * ( v_ice(ji-1,jj) + v_ice(ji,jj) ) 
2413                        END DO
2414                     END DO
2415#endif
2416                  ENDIF
2417               END SELECT
2418            END SELECT
2419            CALL lbc_lnk( zotx1, ssnd(jps_ocx1)%clgrid, -1. )   ;   CALL lbc_lnk( zoty1, ssnd(jps_ocy1)%clgrid, -1. )
2420            !
2421         ENDIF
2422         !
2423         !
2424         IF( TRIM( sn_snd_crt%clvor ) == 'eastward-northward' ) THEN             ! Rotation of the components
2425            !                                                                     ! Ocean component
2426            IF ( TRIM( sn_snd_crt%clvgrd ) == 'T' ) THEN
2427               CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->e', ztmp1 )       ! 1st component
2428               CALL rot_rep( zotx1, zoty1, ssnd(jps_ocx1)%clgrid, 'ij->n', ztmp2 )       ! 2nd component
2429               zotx1(:,:) = ztmp1(:,:)                                                   ! overwrite the components
2430               zoty1(:,:) = ztmp2(:,:)
2431               IF( ssnd(jps_ivx1)%laction ) THEN                                  ! Ice component
2432                  CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->e', ztmp1 )    ! 1st component
2433                  CALL rot_rep( zitx1, zity1, ssnd(jps_ivx1)%clgrid, 'ij->n', ztmp2 )    ! 2nd component
2434                  zitx1(:,:) = ztmp1(:,:)                                                ! overwrite the components
2435                  zity1(:,:) = ztmp2(:,:)
2436               ENDIF
2437            ELSE
2438               ! Temporary code for HadGEM3 - will be removed eventually.
2439               ! Only applies when we want uvel on U grid and vvel on V grid
2440               ! Rotate U and V onto geographic grid before sending.
2441
2442               DO jj=2,jpjm1
2443                  DO ji=2,jpim1
2444                     ztmp1(ji,jj)=0.25*vmask(ji,jj,1)                  &
2445                          *(zotx1(ji,jj)+zotx1(ji-1,jj)    &
2446                          +zotx1(ji,jj+1)+zotx1(ji-1,jj+1))
2447                     ztmp2(ji,jj)=0.25*umask(ji,jj,1)                  &
2448                          *(zoty1(ji,jj)+zoty1(ji+1,jj)    &
2449                          +zoty1(ji,jj-1)+zoty1(ji+1,jj-1))
2450                  ENDDO
2451               ENDDO
2452               
2453               ! Ensure any N fold and wrap columns are updated
2454               CALL lbc_lnk(ztmp1, 'V', -1.0)
2455               CALL lbc_lnk(ztmp2, 'U', -1.0)
2456               
2457               ikchoix = -1
2458               CALL repcmo (zotx1,ztmp2,ztmp1,zoty1,zotx1,zoty1,ikchoix)
2459           ENDIF
2460         ENDIF
2461         !
2462         ! spherical coordinates to cartesian -> 2 components to 3 components
2463         IF( TRIM( sn_snd_crt%clvref ) == 'cartesian' ) THEN
2464            ztmp1(:,:) = zotx1(:,:)                     ! ocean currents
2465            ztmp2(:,:) = zoty1(:,:)
2466            CALL oce2geo ( ztmp1, ztmp2, 'T', zotx1, zoty1, zotz1 )
2467            !
2468            IF( ssnd(jps_ivx1)%laction ) THEN           ! ice velocities
2469               ztmp1(:,:) = zitx1(:,:)
2470               ztmp1(:,:) = zity1(:,:)
2471               CALL oce2geo ( ztmp1, ztmp2, 'T', zitx1, zity1, zitz1 )
2472            ENDIF
2473         ENDIF
2474         !
2475         IF( ssnd(jps_ocx1)%laction )   CALL cpl_snd( jps_ocx1, isec, RESHAPE ( zotx1, (/jpi,jpj,1/) ), info )   ! ocean x current 1st grid
2476         IF( ssnd(jps_ocy1)%laction )   CALL cpl_snd( jps_ocy1, isec, RESHAPE ( zoty1, (/jpi,jpj,1/) ), info )   ! ocean y current 1st grid
2477         IF( ssnd(jps_ocz1)%laction )   CALL cpl_snd( jps_ocz1, isec, RESHAPE ( zotz1, (/jpi,jpj,1/) ), info )   ! ocean z current 1st grid
2478         !
2479         IF( ssnd(jps_ivx1)%laction )   CALL cpl_snd( jps_ivx1, isec, RESHAPE ( zitx1, (/jpi,jpj,1/) ), info )   ! ice   x current 1st grid
2480         IF( ssnd(jps_ivy1)%laction )   CALL cpl_snd( jps_ivy1, isec, RESHAPE ( zity1, (/jpi,jpj,1/) ), info )   ! ice   y current 1st grid
2481         IF( ssnd(jps_ivz1)%laction )   CALL cpl_snd( jps_ivz1, isec, RESHAPE ( zitz1, (/jpi,jpj,1/) ), info )   ! ice   z current 1st grid
2482         !
2483      ENDIF
2484      !
2485      !
2486      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling
2487      !                                                        ! SSH
2488      IF( ssnd(jps_ssh )%laction )  THEN
2489         !                          ! removed inverse barometer ssh when Patm
2490         !                          forcing is used (for sea-ice dynamics)
2491         IF( ln_apr_dyn ) THEN   ;   ztmp1(:,:) = sshb(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
2492         ELSE                    ;   ztmp1(:,:) = sshn(:,:)
2493         ENDIF
2494         CALL cpl_snd( jps_ssh   , isec, RESHAPE ( ztmp1            , (/jpi,jpj,1/) ), info )
2495
2496      ENDIF
2497      !                                                        ! SSS
2498      IF( ssnd(jps_soce  )%laction )  THEN
2499         CALL cpl_snd( jps_soce  , isec, RESHAPE ( tsn(:,:,1,jp_sal), (/jpi,jpj,1/) ), info )
2500      ENDIF
2501      !                                                        ! first T level thickness
2502      IF( ssnd(jps_e3t1st )%laction )  THEN
2503         CALL cpl_snd( jps_e3t1st, isec, RESHAPE ( fse3t_n(:,:,1)   , (/jpi,jpj,1/) ), info )
2504      ENDIF
2505      !                                                        ! Qsr fraction
2506      IF( ssnd(jps_fraqsr)%laction )  THEN
2507         CALL cpl_snd( jps_fraqsr, isec, RESHAPE ( fraqsr_1lev(:,:) , (/jpi,jpj,1/) ), info )
2508      ENDIF
2509      !
2510      !  Fields sent by SAS to OPA when OASIS coupling
2511      !                                                        ! Solar heat flux
2512      IF( ssnd(jps_qsroce)%laction )  CALL cpl_snd( jps_qsroce, isec, RESHAPE ( qsr , (/jpi,jpj,1/) ), info )
2513      IF( ssnd(jps_qnsoce)%laction )  CALL cpl_snd( jps_qnsoce, isec, RESHAPE ( qns , (/jpi,jpj,1/) ), info )
2514      IF( ssnd(jps_oemp  )%laction )  CALL cpl_snd( jps_oemp  , isec, RESHAPE ( emp , (/jpi,jpj,1/) ), info )
2515      IF( ssnd(jps_sflx  )%laction )  CALL cpl_snd( jps_sflx  , isec, RESHAPE ( sfx , (/jpi,jpj,1/) ), info )
2516      IF( ssnd(jps_otx1  )%laction )  CALL cpl_snd( jps_otx1  , isec, RESHAPE ( utau, (/jpi,jpj,1/) ), info )
2517      IF( ssnd(jps_oty1  )%laction )  CALL cpl_snd( jps_oty1  , isec, RESHAPE ( vtau, (/jpi,jpj,1/) ), info )
2518      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info )
2519      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info )
2520     
2521      ztmp1(:,:) = sstfrz(:,:) + rt0
2522      IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
2523      !
2524      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 )
2525      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 )
2526      !
2527      IF( nn_timing == 1 )  CALL timing_stop('sbc_cpl_snd')
2528      !
2529   END SUBROUTINE sbc_cpl_snd
2530   
2531   !!======================================================================
2532END MODULE sbccpl
Note: See TracBrowser for help on using the repository browser.