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

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

Remove RSRH print statements.

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