New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

Last change on this file was 13070, checked in by timgraham, 4 years ago

Added control over printing of iceshelf mass information

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