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

Last change on this file since 12577 was 12577, checked in by dancopsey, 22 months ago

Add a couple of missing changes added through changesets 7540 and 8046 of the original GO6 package branch dev_r5518_GO6_package

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