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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.

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