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

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

sbcice_cice.F90 in branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8198

Last change on this file since 8198 was 8198, checked in by jwhile, 7 years ago
Small fix to sbcice_cice
File size: 47.5 KB

Line
1	MODULE sbcice_cice
2	!!======================================================================
3	!! * MODULE sbcice_cice *
4	!! To couple with sea ice model CICE (LANL)
5	!!=====================================================================
6	#if defined key_cice
7	!!----------------------------------------------------------------------
8	!! 'key_cice' : CICE sea-ice model
9	!!----------------------------------------------------------------------
10	!! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area
11	!!
12	!!
13	!!----------------------------------------------------------------------
14	USE oce ! ocean dynamics and tracers
15	USE dom_oce ! ocean space and time domain
16	USE domvvl
17	USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
18	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
19	USE in_out_manager ! I/O manager
20	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
21	USE lib_mpp ! distributed memory computing library
22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
23	USE wrk_nemo ! work arrays
24	USE timing ! Timing
25	USE daymod ! calendar
26	USE fldread ! read input fields
27	USE sbc_oce ! Surface boundary condition: ocean fields
28	USE sbc_ice ! Surface boundary condition: ice fields
29	USE sbcblk_core ! Surface boundary condition: CORE bulk
30	USE sbccpl
31
32	USE ice_kinds_mod
33	USE ice_blocks
34	USE ice_domain
35	USE ice_domain_size
36	USE ice_boundary
37	USE ice_constants
38	USE ice_gather_scatter
39	USE ice_calendar, only: dt
40	# if defined key_cice4
41	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
42	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
43	strocnxT,strocnyT, &
44	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
45	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
46	flatn_f,fsurfn_f,fcondtopn_f, &
47	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
48	swvdr,swvdf,swidr,swidf,Tf
49	USE ice_therm_vertical, only: calc_Tsfc
50	#else
51	USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
52	vsnon,vice,vicen,nt_Tsfc
53	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
54	strocnxT,strocnyT, &
55	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
56	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
57	flatn_f,fsurfn_f,fcondtopn_f, &
58	#ifdef key_asminc
59	daice_da,fresh_da,fsalt_da, &
60	#endif
61	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
62	swvdr,swvdf,swidr,swidf,Tf, &
63	!! When using NEMO with CICE, this change requires use of
64	!! one of the following two CICE branches:
65	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
66	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
67	keffn_top,Tn_top
68
69	USE ice_therm_shared, only: calc_Tsfc, heat_capacity
70	USE ice_shortwave, only: apeffn
71	#endif
72	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
73	USE ice_atmo, only: calc_strair
74
75	USE CICE_InitMod
76	USE CICE_RunMod
77	USE CICE_FinalMod
78
79	IMPLICIT NONE
80	PRIVATE
81
82	!! * Routine accessibility
83	PUBLIC cice_sbc_init ! routine called by sbc_init
84	PUBLIC cice_sbc_final ! routine called by sbc_final
85	PUBLIC sbc_ice_cice ! routine called by sbc
86
87	INTEGER :: ji_off
88	INTEGER :: jj_off
89
90	INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read
91	INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file
92	INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file
93	INTEGER , PARAMETER :: jp_sblm = 3 ! index of sublimation file
94	INTEGER , PARAMETER :: jp_top1 = 4 ! index of category 1 topmelt file
95	INTEGER , PARAMETER :: jp_top2 = 5 ! index of category 2 topmelt file
96	INTEGER , PARAMETER :: jp_top3 = 6 ! index of category 3 topmelt file
97	INTEGER , PARAMETER :: jp_top4 = 7 ! index of category 4 topmelt file
98	INTEGER , PARAMETER :: jp_top5 = 8 ! index of category 5 topmelt file
99	INTEGER , PARAMETER :: jp_bot1 = 9 ! index of category 1 botmelt file
100	INTEGER , PARAMETER :: jp_bot2 = 10 ! index of category 2 botmelt file
101	INTEGER , PARAMETER :: jp_bot3 = 11 ! index of category 3 botmelt file
102	INTEGER , PARAMETER :: jp_bot4 = 12 ! index of category 4 botmelt file
103	INTEGER , PARAMETER :: jp_bot5 = 13 ! index of category 5 botmelt file
104	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
105
106	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice
107
108	!! * Substitutions
109	# include "domzgr_substitute.h90"
110
111	!! $Id$
112	CONTAINS
113
114	INTEGER FUNCTION sbc_ice_cice_alloc()
115	!!----------------------------------------------------------------------
116	!! * FUNCTION sbc_ice_cice_alloc *
117	!!----------------------------------------------------------------------
118	ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
119	IF( lk_mpp ) CALL mpp_sum ( sbc_ice_cice_alloc )
120	IF( sbc_ice_cice_alloc > 0 ) CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
121	END FUNCTION sbc_ice_cice_alloc
122
123	SUBROUTINE sbc_ice_cice( kt, ksbc )
124	!!---------------------------------------------------------------------
125	!! * ROUTINE sbc_ice_cice *
126	!!
127	!! ** Purpose : update the ocean surface boundary condition via the
128	!! CICE Sea Ice Model time stepping
129	!!
130	!! ** Method : - Get any extra forcing fields for CICE
131	!! - Prepare forcing fields
132	!! - CICE model time stepping
133	!! - call the routine that computes mass and
134	!! heat fluxes at the ice/ocean interface
135	!!
136	!! ** Action : - time evolution of the CICE sea-ice model
137	!! - update all sbc variables below sea-ice:
138	!! utau, vtau, qns , qsr, emp , sfx
139	!!---------------------------------------------------------------------
140	INTEGER, INTENT(in) :: kt ! ocean time step
141	INTEGER, INTENT(in) :: ksbc ! surface forcing type
142	!!----------------------------------------------------------------------
143	!
144	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_cice')
145	!
146	! !----------------------!
147	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
148	! !----------------------!
149
150	! Make sure any fluxes required for CICE are set
151	IF ( ksbc == jp_flx ) THEN
152	CALL cice_sbc_force(kt)
153	ELSE IF ( ksbc == jp_purecpl ) THEN
154	CALL sbc_cpl_ice_flx( 1.0-fr_i )
155	ENDIF
156
157	CALL cice_sbc_in ( kt, ksbc )
158	CALL CICE_Run
159	CALL cice_sbc_out ( kt, ksbc )
160
161	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(kt+1)
162
163	ENDIF ! End sea-ice time step only
164	!
165	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_cice')
166
167	END SUBROUTINE sbc_ice_cice
168
169	SUBROUTINE cice_sbc_init (ksbc)
170	!!---------------------------------------------------------------------
171	!! * ROUTINE cice_sbc_init *
172	!! ** Purpose: Initialise ice related fields for NEMO and coupling
173	!!
174	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
175	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
176	REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d
177	INTEGER :: ji, jj, jl, jk ! dummy loop indices
178	!!---------------------------------------------------------------------
179
180	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_init')
181	!
182	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
183	!
184	IF(lwp) WRITE(numout,*)'cice_sbc_init'
185
186	ji_off = INT ( (jpiglo - nx_global) / 2 )
187	jj_off = INT ( (jpjglo - ny_global) / 2 )
188
189	! Initialize CICE
190	CALL CICE_Initialize
191
192	! Do some CICE consistency checks
193	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
194	IF ( calc_strair .OR. calc_Tsfc ) THEN
195	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
196	ENDIF
197	ELSEIF (ksbc == jp_core) THEN
198	IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
199	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
200	ENDIF
201	ENDIF
202
203
204	! allocate sbc_ice and sbc_cice arrays
205	IF( sbc_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' )
206	IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
207
208	! Ensure that no temperature points are below freezing if not a NEMO restart
209	IF( .NOT. ln_rstart ) THEN
210
211	CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d )
212	DO jk=1,jpk
213	CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), fsdept_n(:,:,jk) )
214	ENDDO
215	tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
216	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
217	CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d )
218
219	#if defined key_nemocice_decomp
220	! Pass initial SST from NEMO to CICE so ice is initialised correctly if
221	! there is no restart file.
222	! Values from a CICE restart file would overwrite this
223	CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)
224	#endif
225
226	ENDIF
227
228	! calculate surface freezing temperature and send to CICE
229	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1))
230	CALL nemo2cice(sstfrz,Tf, 'T', 1. )
231
232	CALL cice2nemo(aice,fr_i, 'T', 1. )
233	IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
234	DO jl=1,ncat
235	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
236	ENDDO
237	ENDIF
238
239	! T point to U point
240	! T point to V point
241	fr_iu(:,:)=0.0
242	fr_iv(:,:)=0.0
243	DO jj=1,jpjm1
244	DO ji=1,jpim1
245	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
246	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
247	ENDDO
248	ENDDO
249
250	CALL lbc_lnk ( fr_iu , 'U', 1. )
251	CALL lbc_lnk ( fr_iv , 'V', 1. )
252
253	! ! embedded sea ice
254	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
255	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
256	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
257	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
258	snwice_mass_b(:,:) = snwice_mass(:,:)
259	ELSE
260	snwice_mass (:,:) = 0.0_wp ! no mass exchanges
261	snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges
262	ENDIF
263	IF( .NOT. ln_rstart ) THEN
264	IF( nn_ice_embd == 2 ) THEN ! full embedment (case 2) deplete the initial ssh below sea-ice area
265	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
266	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
267	#if defined key_vvl
268	! key_vvl necessary? clem: yes for compilation purpose
269	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
270	fse3t_n(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshn(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
271	fse3t_b(:,:,jk) = e3t_0(:,:,jk)( 1._wp + sshb(:,:)tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
272	ENDDO
273	fse3t_a(:,:,:) = fse3t_b(:,:,:)
274	! Reconstruction of all vertical scale factors at now and before time
275	! steps
276	! =============================================================================
277	! Horizontal scale factor interpolations
278	! --------------------------------------
279	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
280	CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
281	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
282	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
283	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
284	! Vertical scale factor interpolations
285	! ------------------------------------
286	CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
287	CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
288	CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
289	CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
290	CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
291	! t- and w- points depth
292	! ----------------------
293	fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
294	fsdepw_n(:,:,1) = 0.0_wp
295	fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
296	DO jk = 2, jpk
297	fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
298	fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
299	fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:)
300	END DO
301	#endif
302	ENDIF
303	ENDIF
304
305	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
306
307	#if defined key_asminc
308	! Initialize fresh water and salt fluxes from data assim
309	! and data assimilation index to cice
310	nfresh_da(:,:) = 0.0
311	nfsalt_da(:,:) = 0.0
312	ndaice_da(:,:) = 0.0
313	#endif
314	!
315	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_init')
316	!
317	END SUBROUTINE cice_sbc_init
318
319
320	SUBROUTINE cice_sbc_in (kt, ksbc)
321	!!---------------------------------------------------------------------
322	!! * ROUTINE cice_sbc_in *
323	!! ** Purpose: Set coupling fields and pass to CICE
324	!!---------------------------------------------------------------------
325	INTEGER, INTENT(in ) :: kt ! ocean time step
326	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
327
328	INTEGER :: ji, jj, jl ! dummy loop indices
329	REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
330	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
331	REAL(wp) :: zintb, zintn ! dummy argument
332	!!---------------------------------------------------------------------
333
334	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_in')
335	!
336	CALL wrk_alloc( jpi,jpj, ztmp, zpice )
337	CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
338
339	IF( kt == nit000 ) THEN
340	IF(lwp) WRITE(numout,*)'cice_sbc_in'
341	ENDIF
342
343	ztmp(:,:)=0.0
344
345	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
346	! the first time-step)
347
348	! forced and coupled case
349
350	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
351
352	ztmpn(:,:,:)=0.0
353
354	! x comp of wind stress (CI_1)
355	! U point to F point
356	DO jj=1,jpjm1
357	DO ji=1,jpi
358	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
359	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
360	ENDDO
361	ENDDO
362	CALL nemo2cice(ztmp,strax,'F', -1. )
363
364	! y comp of wind stress (CI_2)
365	! V point to F point
366	DO jj=1,jpj
367	DO ji=1,jpim1
368	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
369	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
370	ENDDO
371	ENDDO
372	CALL nemo2cice(ztmp,stray,'F', -1. )
373
374
375	! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
376	! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby
377	! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
378	! gridbox mean fluxes in the UM by future ice concentration obtained through
379	! OASIS. This allows for a much more realistic apportionment of energy through
380	! the ice - and conserves energy.
381	! Therefore the fluxes are now divided by ice concentration in the coupled
382	! formulation (jp_purecpl) as well as for jp_flx. This NEMO branch should only
383	! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
384	! which point the GSI8 UM changes were committed.
385
386	! Surface downward latent heat flux (CI_5)
387	IF (ksbc == jp_flx) THEN
388	DO jl=1,ncat
389	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
390	ENDDO
391	ELSE IF (ksbc == jp_purecpl) THEN
392	DO jl=1,ncat
393	ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
394	ENDDO
395	ELSE
396	!In coupled mode - qla_ice calculated in sbc_cpl for each category
397	ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
398	ENDIF
399
400	DO jl=1,ncat
401	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
402
403	! GBM conductive flux through ice (CI_6)
404	! Convert to GBM
405	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
406	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
407	ELSE
408	ztmp(:,:) = botmelt(:,:,jl)
409	ENDIF
410	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
411
412	! GBM surface heat flux (CI_7)
413	! Convert to GBM
414	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
415	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
416	ELSE
417	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
418	ENDIF
419	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
420	ENDDO
421
422	ELSE IF (ksbc == jp_core) THEN
423
424	! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
425	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
426	ztmp(:,:) = wndi_ice(:,:)
427	CALL nemo2cice(ztmp,uatm,'T', -1. )
428	ztmp(:,:) = wndj_ice(:,:)
429	CALL nemo2cice(ztmp,vatm,'T', -1. )
430	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
431	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
432	ztmp(:,:) = qsr_ice(:,:,1)
433	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
434	ztmp(:,:) = qlw_ice(:,:,1)
435	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
436	ztmp(:,:) = tatm_ice(:,:)
437	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
438	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
439	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
440	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
441	! Constant (101000.) atm pressure assumed
442	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
443	ztmp(:,:) = qatm_ice(:,:)
444	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
445	ztmp(:,:)=10.0
446	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
447
448	! May want to check all values are physically realistic (as in CICE routine
449	! prepare_forcing)?
450
451	! Divide shortwave into spectral bands (as in prepare_forcing)
452	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
453	CALL nemo2cice(ztmp,swvdr,'T', 1. )
454	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
455	CALL nemo2cice(ztmp,swvdf,'T', 1. )
456	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
457	CALL nemo2cice(ztmp,swidr,'T', 1. )
458	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
459	CALL nemo2cice(ztmp,swidf,'T', 1. )
460
461	ENDIF
462
463	#if defined key_asminc
464	!Ice concentration change (from assimilation)
465	ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
466	Call nemo2cice(ztmp,daice_da,'T', 1. )
467	#endif
468
469	! Snowfall
470	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
471	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
472	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
473	CALL nemo2cice(ztmp,fsnow,'T', 1. )
474
475	! Rainfall
476	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
477	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
478	CALL nemo2cice(ztmp,frain,'T', 1. )
479
480	! Recalculate freezing temperature and send to CICE
481	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1))
482	CALL nemo2cice(sstfrz,Tf,'T', 1. )
483
484	! Freezing/melting potential
485	! Calculated over NEMO leapfrog timestep (hence 2*dt)
486	nfrzmlt(:,:)=rau0rcpfse3t_m(:,:)(sstfrz(:,:)-sst_m(:,:))/(2.0dt)
487	CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
488
489	! SST and SSS
490
491	CALL nemo2cice(sst_m,sst,'T', 1. )
492	CALL nemo2cice(sss_m,sss,'T', 1. )
493
494	IF( ksbc == jp_purecpl ) THEN
495	! Sea ice surface skin temperature
496	DO jl=1,ncat
497	CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
498	ENDDO
499	ENDIF
500
501	! x comp and y comp of surface ocean current
502	! U point to F point
503	DO jj=1,jpjm1
504	DO ji=1,jpi
505	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
506	ENDDO
507	ENDDO
508	CALL nemo2cice(ztmp,uocn,'F', -1. )
509
510	! V point to F point
511	DO jj=1,jpj
512	DO ji=1,jpim1
513	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
514	ENDDO
515	ENDDO
516	CALL nemo2cice(ztmp,vocn,'F', -1. )
517
518	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
519	!
520	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
521	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
522	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
523	!
524	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
525	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
526	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
527	!
528	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
529	!
530	!
531	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
532	zpice(:,:) = ssh_m(:,:)
533	ENDIF
534
535	! x comp and y comp of sea surface slope (on F points)
536	! T point to F point
537	DO jj=1,jpjm1
538	DO ji=1,jpim1
539	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj ))/e1u(ji,jj ) &
540	+ (zpice(ji+1,jj+1)-zpice(ji,jj+1))/e1u(ji,jj+1) ) &
541	* fmask(ji,jj,1)
542	ENDDO
543	ENDDO
544	CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
545
546	! T point to F point
547	DO jj=1,jpjm1
548	DO ji=1,jpim1
549	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj))/e2v(ji ,jj) &
550	+ (zpice(ji+1,jj+1)-zpice(ji+1,jj))/e2v(ji+1,jj) ) &
551	* fmask(ji,jj,1)
552	ENDDO
553	ENDDO
554	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
555
556	CALL wrk_dealloc( jpi,jpj, ztmp, zpice )
557	CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
558	!
559	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_in')
560	!
561	END SUBROUTINE cice_sbc_in
562
563
564	SUBROUTINE cice_sbc_out (kt,ksbc)
565	!!---------------------------------------------------------------------
566	!! * ROUTINE cice_sbc_out *
567	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
568	!!---------------------------------------------------------------------
569	INTEGER, INTENT( in ) :: kt ! ocean time step
570	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
571
572	INTEGER :: ji, jj, jl ! dummy loop indices
573	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
574	!!---------------------------------------------------------------------
575
576	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_out')
577	!
578	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
579
580	IF( kt == nit000 ) THEN
581	IF(lwp) WRITE(numout,*)'cice_sbc_out'
582	ENDIF
583
584	! x comp of ocean-ice stress
585	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
586	ss_iou(:,:)=0.0
587	! F point to U point
588	DO jj=2,jpjm1
589	DO ji=2,jpim1
590	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
591	ENDDO
592	ENDDO
593	CALL lbc_lnk( ss_iou , 'U', -1. )
594
595	! y comp of ocean-ice stress
596	CALL cice2nemo(strocny,ztmp1,'F', -1. )
597	ss_iov(:,:)=0.0
598	! F point to V point
599
600	DO jj=1,jpjm1
601	DO ji=2,jpim1
602	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
603	ENDDO
604	ENDDO
605	CALL lbc_lnk( ss_iov , 'V', -1. )
606
607	! x and y comps of surface stress
608	! Combine wind stress and ocean-ice stress
609	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
610	! strocnx and strocny already weighted by ice fraction in CICE so not done here
611
612	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
613	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
614
615	! Also need ice/ocean stress on T points so that taum can be updated
616	! This interpolation is already done in CICE so best to use those values
617	CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
618	CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
619
620	! Update taum with modulus of ice-ocean stress
621	! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
622	taum(:,:)=(1.0-fr_i(:,:))taum(:,:)+fr_i(:,:)SQRT(ztmp12. + ztmp22.)
623
624	! Freshwater fluxes
625
626	IF (ksbc == jp_flx) THEN
627	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
628	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
629	! Not ideal since aice won't be the same as in the atmosphere.
630	! Better to use evap and tprecip? (but for now don't read in evap in this case)
631	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
632	ELSE IF (ksbc == jp_core) THEN
633	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
634	ELSE IF (ksbc == jp_purecpl) THEN
635	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
636	! This is currently as required with the coupling fields from the UM atmosphere
637	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
638	ENDIF
639
640	#if defined key_cice4
641	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
642	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
643	#else
644	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
645	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
646	#endif
647
648	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
649	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
650	! which has to be compensated for by the ocean salinity potentially going negative
651	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
652	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
653	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
654	sfx(:,:)=ztmp2(:,:)*1000.0
655	emp(:,:)=emp(:,:)-ztmp1(:,:)
656	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
657
658	CALL lbc_lnk( emp , 'T', 1. )
659	CALL lbc_lnk( sfx , 'T', 1. )
660
661	! Solar penetrative radiation and non solar surface heat flux
662
663	! Scale qsr and qns according to ice fraction (bulk formulae only)
664
665	IF (ksbc == jp_core) THEN
666	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
667	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
668	ENDIF
669	! Take into account snow melting except for fully coupled when already in qns_tot
670	IF (ksbc == jp_purecpl) THEN
671	qsr(:,:)= qsr_tot(:,:)
672	qns(:,:)= qns_tot(:,:)
673	ELSE
674	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
675	ENDIF
676
677	! Now add in ice / snow related terms
678	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
679	#if defined key_cice4
680	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
681	#else
682	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
683	#endif
684	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
685	CALL lbc_lnk( qsr , 'T', 1. )
686
687	DO jj=1,jpj
688	DO ji=1,jpi
689	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
690	ENDDO
691	ENDDO
692
693	#if defined key_cice4
694	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
695	#else
696	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
697	#endif
698	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
699
700	CALL lbc_lnk( qns , 'T', 1. )
701
702	! Prepare for the following CICE time-step
703
704	CALL cice2nemo(aice,fr_i,'T', 1. )
705	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
706	DO jl=1,ncat
707	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
708	ENDDO
709	ENDIF
710
711	! T point to U point
712	! T point to V point
713	DO jj=1,jpjm1
714	DO ji=1,jpim1
715	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
716	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
717	ENDDO
718	ENDDO
719
720	CALL lbc_lnk ( fr_iu , 'U', 1. )
721	CALL lbc_lnk ( fr_iv , 'V', 1. )
722
723	! ! embedded sea ice
724	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
725	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
726	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
727	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
728	snwice_mass_b(:,:) = snwice_mass(:,:)
729	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
730	ENDIF
731
732	#if defined key_asminc
733	! Import fresh water and salt flux due to seaice da
734	CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
735	CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
736	#endif
737
738	! Release work space
739
740	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
741	!
742	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_out')
743	!
744	END SUBROUTINE cice_sbc_out
745
746
747	SUBROUTINE cice_sbc_hadgam( kt )
748	!!---------------------------------------------------------------------
749	!! * ROUTINE cice_sbc_hadgam *
750	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
751	!!
752	!!
753	INTEGER, INTENT( in ) :: kt ! ocean time step
754	!!---------------------------------------------------------------------
755
756	INTEGER :: jl ! dummy loop index
757	INTEGER :: ierror
758
759	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_hadgam')
760	!
761	IF( kt == nit000 ) THEN
762	IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
763	IF( sbc_cpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
764	ENDIF
765
766	! ! =========================== !
767	! ! Prepare Coupling fields !
768	! ! =========================== !
769
770	! x and y comp of ice velocity
771
772	CALL cice2nemo(uvel,u_ice,'F', -1. )
773	CALL cice2nemo(vvel,v_ice,'F', -1. )
774
775	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
776
777	! Snow and ice thicknesses (CO_2 and CO_3)
778
779	DO jl = 1,ncat
780	CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
781	CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
782	ENDDO
783
784	#if ! defined key_cice4
785	! Meltpond fraction and depth
786	DO jl = 1,ncat
787	CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. )
788	CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. )
789	ENDDO
790	#endif
791
792
793	! If using multilayers thermodynamics in CICE then get top layer temperature
794	! and effective conductivity
795	!! When using NEMO with CICE, this change requires use of
796	!! one of the following two CICE branches:
797	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
798	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
799	IF (heat_capacity) THEN
800	DO jl = 1,ncat
801	CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. )
802	CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. )
803	ENDDO
804	! Convert surface temperature to Kelvin
805	tn_ice(:,:,:)=tn_ice(:,:,:)+rt0
806	ELSE
807	tn_ice(:,:,:) = 0.0
808	kn_ice(:,:,:) = 0.0
809	ENDIF
810
811	!
812	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
813	!
814	END SUBROUTINE cice_sbc_hadgam
815
816
817	SUBROUTINE cice_sbc_final
818	!!---------------------------------------------------------------------
819	!! * ROUTINE cice_sbc_final *
820	!! ** Purpose: Finalize CICE
821	!!---------------------------------------------------------------------
822
823	IF(lwp) WRITE(numout,*)'cice_sbc_final'
824
825	CALL CICE_Finalize
826
827	END SUBROUTINE cice_sbc_final
828
829	SUBROUTINE cice_sbc_force (kt)
830	!!---------------------------------------------------------------------
831	!! * ROUTINE cice_sbc_force *
832	!! ** Purpose : Provide CICE forcing from files
833	!!
834	!!---------------------------------------------------------------------
835	!! ** Method : READ monthly flux file in NetCDF files
836	!!
837	!! snowfall
838	!! rainfall
839	!! sublimation rate
840	!! topmelt (category)
841	!! botmelt (category)
842	!!
843	!! History :
844	!!----------------------------------------------------------------------
845	!! * Modules used
846	USE iom
847
848	!! * arguments
849	INTEGER, INTENT( in ) :: kt ! ocean time step
850
851	INTEGER :: ierror ! return error code
852	INTEGER :: ifpr ! dummy loop index
853	!!
854	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
855	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
856	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
857	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
858	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
859
860	!!
861	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
862	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
863	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
864	INTEGER :: ios
865	!!---------------------------------------------------------------------
866
867	! ! ====================== !
868	IF( kt == nit000 ) THEN ! First call kt=nit000 !
869	! ! ====================== !
870	! namsbc_cice is not yet in the reference namelist
871	! set file information (default values)
872	cn_dir = './' ! directory in which the model is executed
873
874	! (NB: frequency positive => hours, negative => months)
875	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
876	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
877	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
878	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
879	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
880	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
881	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
882	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
883	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
884	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
885	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
886	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
887	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
888	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
889	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
890
891	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
892	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
893	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
894
895	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
896	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
897	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
898	IF(lwm) WRITE ( numond, namsbc_cice )
899
900	! store namelist information in an array
901	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
902	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
903	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
904	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
905	slf_i(jp_bot5) = sn_bot5
906
907	! set sf structure
908	ALLOCATE( sf(jpfld), STAT=ierror )
909	IF( ierror > 0 ) THEN
910	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
911	ENDIF
912
913	DO ifpr= 1, jpfld
914	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
915	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
916	END DO
917
918	! fill sf with slf_i and control print
919	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
920	!
921	ENDIF
922
923	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
924	! ! input fields at the current time-step
925
926	! set the fluxes from read fields
927	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
928	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
929	! May be better to do this conversion somewhere else
930	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
931	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
932	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
933	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
934	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
935	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
936	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
937	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
938	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
939	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
940	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
941
942	! control print (if less than 100 time-step asked)
943	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
944	WRITE(numout,*)
945	WRITE(numout,*) ' read forcing fluxes for CICE OK'
946	CALL FLUSH(numout)
947	ENDIF
948
949	END SUBROUTINE cice_sbc_force
950
951	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
952	!!---------------------------------------------------------------------
953	!! * ROUTINE nemo2cice *
954	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
955	#if defined key_nemocice_decomp
956	!!
957	!! NEMO and CICE PE sub domains are identical, hence
958	!! there is no need to gather or scatter data from
959	!! one PE configuration to another.
960	#else
961	!! Automatically gather/scatter between
962	!! different processors and blocks
963	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
964	!! B. Gather pn into global array (png)
965	!! C. Map png into CICE global array (pcg)
966	!! D. Scatter pcg to CICE blocks (pc) + update haloes
967	#endif
968	!!---------------------------------------------------------------------
969
970	CHARACTER(len=1), INTENT( in ) :: &
971	cd_type ! nature of pn grid-point
972	! ! = T or F gridpoints
973	REAL(wp), INTENT( in ) :: &
974	psgn ! control of the sign change
975	! ! =-1 , the sign is modified following the type of b.c. used
976	! ! = 1 , no sign change
977	REAL(wp), DIMENSION(jpi,jpj) :: pn
978	#if !defined key_nemocice_decomp
979	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
980	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
981	#endif
982	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
983	INTEGER (int_kind) :: &
984	field_type, &! id for type of field (scalar, vector, angle)
985	grid_loc ! id for location on horizontal grid
986	! (center, NEcorner, Nface, Eface)
987
988	INTEGER :: ji, jj, jn ! dummy loop indices
989
990	! A. Ensure all haloes are filled in NEMO field (pn)
991
992	CALL lbc_lnk( pn , cd_type, psgn )
993
994	#if defined key_nemocice_decomp
995
996	! Copy local domain data from NEMO to CICE field
997	pc(:,:,1)=0.0
998	DO jj=2,ny_block-1
999	DO ji=2,nx_block-1
1000	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
1001	ENDDO
1002	ENDDO
1003
1004	#else
1005
1006	! B. Gather pn into global array (png)
1007
1008	IF ( jpnij > 1) THEN
1009	CALL mppsync
1010	CALL mppgather (pn,0,png)
1011	CALL mppsync
1012	ELSE
1013	png(:,:,1)=pn(:,:)
1014	ENDIF
1015
1016	! C. Map png into CICE global array (pcg)
1017
1018	! Need to make sure this is robust to changes in NEMO halo rows....
1019	! (may be OK but not 100% sure)
1020
1021	IF (nproc==0) THEN
1022	! pcg(:,:)=0.0
1023	DO jn=1,jpnij
1024	DO jj=nldjt(jn),nlejt(jn)
1025	DO ji=nldit(jn),nleit(jn)
1026	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
1027	ENDDO
1028	ENDDO
1029	ENDDO
1030	DO jj=1,ny_global
1031	DO ji=1,nx_global
1032	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
1033	ENDDO
1034	ENDDO
1035	ENDIF
1036
1037	#endif
1038
1039	SELECT CASE ( cd_type )
1040	CASE ( 'T' )
1041	grid_loc=field_loc_center
1042	CASE ( 'F' )
1043	grid_loc=field_loc_NEcorner
1044	END SELECT
1045
1046	SELECT CASE ( NINT(psgn) )
1047	CASE ( -1 )
1048	field_type=field_type_vector
1049	CASE ( 1 )
1050	field_type=field_type_scalar
1051	END SELECT
1052
1053	#if defined key_nemocice_decomp
1054	! Ensure CICE halos are up to date
1055	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1056	#else
1057	! D. Scatter pcg to CICE blocks (pc) + update halos
1058	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
1059	#endif
1060
1061	END SUBROUTINE nemo2cice
1062
1063	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
1064	!!---------------------------------------------------------------------
1065	!! * ROUTINE cice2nemo *
1066	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
1067	#if defined key_nemocice_decomp
1068	!!
1069	!! NEMO and CICE PE sub domains are identical, hence
1070	!! there is no need to gather or scatter data from
1071	!! one PE configuration to another.
1072	#else
1073	!! Automatically deal with scatter/gather between
1074	!! different processors and blocks
1075	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
1076	!! B. Map pcg into NEMO global array (png)
1077	!! C. Scatter png into NEMO field (pn) for each processor
1078	!! D. Ensure all haloes are filled in pn
1079	#endif
1080	!!---------------------------------------------------------------------
1081
1082	CHARACTER(len=1), INTENT( in ) :: &
1083	cd_type ! nature of pn grid-point
1084	! ! = T or F gridpoints
1085	REAL(wp), INTENT( in ) :: &
1086	psgn ! control of the sign change
1087	! ! =-1 , the sign is modified following the type of b.c. used
1088	! ! = 1 , no sign change
1089	REAL(wp), DIMENSION(jpi,jpj) :: pn
1090
1091	#if defined key_nemocice_decomp
1092	INTEGER (int_kind) :: &
1093	field_type, & ! id for type of field (scalar, vector, angle)
1094	grid_loc ! id for location on horizontal grid
1095	! (center, NEcorner, Nface, Eface)
1096	#else
1097	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
1098	#endif
1099
1100	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
1101
1102	INTEGER :: ji, jj, jn ! dummy loop indices
1103
1104
1105	#if defined key_nemocice_decomp
1106
1107	SELECT CASE ( cd_type )
1108	CASE ( 'T' )
1109	grid_loc=field_loc_center
1110	CASE ( 'F' )
1111	grid_loc=field_loc_NEcorner
1112	END SELECT
1113
1114	SELECT CASE ( NINT(psgn) )
1115	CASE ( -1 )
1116	field_type=field_type_vector
1117	CASE ( 1 )
1118	field_type=field_type_scalar
1119	END SELECT
1120
1121	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1122
1123
1124	pn(:,:)=0.0
1125	DO jj=1,jpjm1
1126	DO ji=1,jpim1
1127	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
1128	ENDDO
1129	ENDDO
1130
1131	#else
1132
1133	! A. Gather CICE blocks (pc) into global array (pcg)
1134
1135	CALL gather_global(pcg, pc, 0, distrb_info)
1136
1137	! B. Map pcg into NEMO global array (png)
1138
1139	! Need to make sure this is robust to changes in NEMO halo rows....
1140	! (may be OK but not spent much time thinking about it)
1141	! Note that non-existent pcg elements may be used below, but
1142	! the lbclnk call on pn will replace these with sensible values
1143
1144	IF (nproc==0) THEN
1145	png(:,:,:)=0.0
1146	DO jn=1,jpnij
1147	DO jj=nldjt(jn),nlejt(jn)
1148	DO ji=nldit(jn),nleit(jn)
1149	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
1150	ENDDO
1151	ENDDO
1152	ENDDO
1153	ENDIF
1154
1155	! C. Scatter png into NEMO field (pn) for each processor
1156
1157	IF ( jpnij > 1) THEN
1158	CALL mppsync
1159	CALL mppscatter (png,0,pn)
1160	CALL mppsync
1161	ELSE
1162	pn(:,:)=png(:,:,1)
1163	ENDIF
1164
1165	#endif
1166
1167	! D. Ensure all haloes are filled in pn
1168
1169	CALL lbc_lnk( pn , cd_type, psgn )
1170
1171	END SUBROUTINE cice2nemo
1172
1173	#else
1174	!!----------------------------------------------------------------------
1175	!! Default option Dummy module NO CICE sea-ice model
1176	!!----------------------------------------------------------------------
1177	!! $Id$
1178	CONTAINS
1179
1180	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
1181	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
1182	END SUBROUTINE sbc_ice_cice
1183
1184	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
1185	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
1186	END SUBROUTINE cice_sbc_init
1187
1188	SUBROUTINE cice_sbc_final ! Dummy routine
1189	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
1190	END SUBROUTINE cice_sbc_final
1191
1192	#endif
1193
1194	!!======================================================================
1195	END MODULE sbcice_cice

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

Download in other formats: