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

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

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10764

Last change on this file since 10764 was 10764, checked in by jcastill, 5 years ago

Minimal set of changes for coupling order

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