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/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8637

Last change on this file since 8637 was 8637, checked in by gm, 7 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626

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