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 @ 8478

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

Moved lines setting up ssnd(jps_fmdice), as they were in the wrong place.

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