source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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