source: branches/UKMO/dev_r5518_GO6_fix_zemp_ice_10681/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 11054

Last change on this file since 11054 was 11054, checked in by dancopsey, 17 months ago

Put the switch back in again.

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