source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10759

Last change on this file since 10759 was 10759, checked in by andmirek, 20 months ago

GMED 450 write output.namelist.dyn only for nprint > 2

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