source: branches/UKMO/dev_r5518_cleanup_1d_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10041

Last change on this file since 10041 was 10041, checked in by dancopsey, 2 years ago

Added coupling of 1D river outflow.

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