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

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

Fix compile issues for when I am running GO6+MEDUSA.

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