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

source: branches/UKMO/dev_r5518_couple_chlorophyll/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8000

Last change on this file since 8000 was 8000, checked in by frrh, 7 years ago

Cater for chlrophyll as a coupling field in the OASIS3-MCT interface.
This data should be populated by MEDUSA, which needs to fill in the
2D T grid point chloro_out_cpl field at the appropriate timesteps.

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