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_r8638_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package_r8638_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 9714

Last change on this file since 9714 was 9714, checked in by jamrae, 6 years ago

Made code changes for sea ice form drag coupling.

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