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_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

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

Last change on this file since 8459 was 8459, checked in by jamrae, 7 years ago

Outputs for debugging.

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