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 @ 8235

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

Update branch in line with current head of GO6 package branch at revision 8104.

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