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

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

source: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 7152

Last change on this file since 7152 was 7152, checked in by jcastill, 7 years ago

Initial implementation of wave coupling branch - INGV wave branch + UKMO wave coupling branch

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