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

source: branches/UKMO/dev_r5518_GO6_fix_zemp_ice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10920

Last change on this file since 10920 was 10920, checked in by dancopsey, 5 years ago

Use ice fractions from the last coupling time to calculate grid box means of coupled fluxes.

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