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

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

source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 8670

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

Added capability to weight sea ice roughness length before coupling.

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