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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbccpl.F90 in branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 10775

Last change on this file since 10775 was 10775, checked in by andmirek, 5 years ago
changes in icesheet/iceberg (coupled configuration)
File size: 157.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: