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

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago

GMED 450 add flush after prints

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