New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbccpl.F90 in branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8385

Last change on this file since 8385 was 8385, checked in by jamrae, 7 years ago

Debugging.

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