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

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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8563

Last change on this file since 8563 was 8563, checked in by clem, 7 years ago

change variable names (ht_s => h_s & ht_i => h_i)

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