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

Last change on this file since 8427 was 8427, checked in by timgraham, 3 years ago

GMED ticket 341. Fixes found through use of preset Na Ns? to make the code more robust.

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