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

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

source: NEMO/branches/UKMO/NEMO_4.0.1_icesheet_and_river_coupling/src/OCE/SBC/sbccpl.F90 @ 12576

Last change on this file since 12576 was 12576, checked in by dancopsey, 4 years ago

Merge in iceberg calving stuff from dev_r5518_coupling_GSI7_GSI8_landice from the start of the branch to revision 6023.

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