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

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

source: branches/UKMO/dev_r5518_GO6_fix_zemp_ice_10681/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 11046

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

Take the switch back out again. It was decided not that this is a bug so should be automatically picked up.

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