source: branches/UKMO/dev_r5518_GO6_package_r7750_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8350

Last change on this file since 8350 was 8350, checked in by jamrae, 4 years ago

Added code to pass sea ice roughness length from CICE to the coupler.

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