Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbccpl.F90 @ 14627

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

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

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

Context Navigation

source: branches/UKMO/jonnywilliams-u-br854/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 @ 14627

Download in other formats: