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_add_pond_lids_prints/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbccpl.F90 @ 12394

Last change on this file since 12394 was 12394, checked in by dancopsey, 5 years ago

Fix compile and XIOS errors.

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