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

Last change on this file since 10775 was 10775, checked in by andmirek, 20 months ago

changes in icesheet/iceberg (coupled configuration)

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