source: branches/UKMO/dev_r5518_GO6_package_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 11274

Last change on this file since 11274 was 11274, checked in by dancopsey, 15 months ago

Merge in sublimation bug fixes (GMED#449).

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