source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8819

Last change on this file since 8819 was 8819, checked in by jamrae, 3 years ago

Implemented code to pass sea ice skin roughness length to atmosphere.

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