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

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

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10478

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

Merged branch UKMO/r6232_HZG_WAVE-coupling@9868

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