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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11883

Last change on this file since 11883 was 11883, checked in by timgraham, 4 years ago
Mergerd in changes from Chris' working copy (except for symlink changes)
File size: 48.2 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	! In coupled mode get extra fields from CICE for passing back to atmosphere
316
317	IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)
318	!
319	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_init')
320	!
321	END SUBROUTINE cice_sbc_init
322
323
324	SUBROUTINE cice_sbc_in (kt, ksbc)
325	!!---------------------------------------------------------------------
326	!! * ROUTINE cice_sbc_in *
327	!! ** Purpose: Set coupling fields and pass to CICE
328	!!---------------------------------------------------------------------
329	INTEGER, INTENT(in ) :: kt ! ocean time step
330	INTEGER, INTENT(in ) :: ksbc ! surface forcing type
331
332	INTEGER :: ji, jj, jl ! dummy loop indices
333	REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
334	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
335	REAL(wp) :: zintb, zintn ! dummy argument
336	!!---------------------------------------------------------------------
337
338	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_in')
339	!
340	CALL wrk_alloc( jpi,jpj, ztmp, zpice )
341	CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
342
343	IF( kt == nit000 ) THEN
344	IF(lwp) WRITE(numout,*)'cice_sbc_in'
345	ENDIF
346
347	ztmp(:,:)=0.0
348
349	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
350	! the first time-step)
351
352	! forced and coupled case
353
354	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
355
356	ztmpn(:,:,:)=0.0
357
358	! x comp of wind stress (CI_1)
359	! U point to F point
360	DO jj=1,jpjm1
361	DO ji=1,jpi
362	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
363	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
364	ENDDO
365	ENDDO
366	CALL nemo2cice(ztmp,strax,'F', -1. )
367
368	! y comp of wind stress (CI_2)
369	! V point to F point
370	DO jj=1,jpj
371	DO ji=1,jpim1
372	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
373	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
374	ENDDO
375	ENDDO
376	CALL nemo2cice(ztmp,stray,'F', -1. )
377
378
379	! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
380	! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby
381	! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
382	! gridbox mean fluxes in the UM by future ice concentration obtained through
383	! OASIS. This allows for a much more realistic apportionment of energy through
384	! the ice - and conserves energy.
385	! Therefore the fluxes are now divided by ice concentration in the coupled
386	! formulation (jp_purecpl) as well as for jp_flx. This NEMO branch should only
387	! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
388	! which point the GSI8 UM changes were committed.
389
390	! Surface downward latent heat flux (CI_5)
391	IF (ksbc == jp_flx) THEN
392	DO jl=1,ncat
393	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
394	ENDDO
395	ELSE IF (ksbc == jp_purecpl) THEN
396	DO jl=1,ncat
397	IF ( ln_fix_sea_ice_fluxes ) THEN
398	ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i_last_couple(:,:,jl)
399	ELSE
400	ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
401	ENDIF
402	ENDDO
403	ELSE
404	!In coupled mode - qla_ice calculated in sbc_cpl for each category
405	ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
406	ENDIF
407
408	DO jl=1,ncat
409	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
410
411	! GBM conductive flux through ice (CI_6)
412	! Convert to GBM
413	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
414	IF ( ln_fix_sea_ice_fluxes ) THEN
415	ztmp(:,:) = botmelt(:,:,jl)*a_i_last_couple(:,:,jl)
416	ELSE
417	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
418	ENDIF
419	ELSE
420	ztmp(:,:) = botmelt(:,:,jl)
421	ENDIF
422	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
423
424	! GBM surface heat flux (CI_7)
425	! Convert to GBM
426	IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
427	IF ( ln_fix_sea_ice_fluxes ) THEN
428	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i_last_couple(:,:,jl)
429	ELSE
430	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
431	ENDIF
432	ELSE
433	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
434	ENDIF
435	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
436	ENDDO
437
438	ELSE IF (ksbc == jp_core) THEN
439
440	! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
441	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
442	ztmp(:,:) = wndi_ice(:,:)
443	CALL nemo2cice(ztmp,uatm,'T', -1. )
444	ztmp(:,:) = wndj_ice(:,:)
445	CALL nemo2cice(ztmp,vatm,'T', -1. )
446	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
447	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
448	ztmp(:,:) = qsr_ice(:,:,1)
449	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
450	ztmp(:,:) = qlw_ice(:,:,1)
451	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
452	ztmp(:,:) = tatm_ice(:,:)
453	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
454	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
455	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
456	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
457	! Constant (101000.) atm pressure assumed
458	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
459	ztmp(:,:) = qatm_ice(:,:)
460	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
461	ztmp(:,:)=10.0
462	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
463
464	! May want to check all values are physically realistic (as in CICE routine
465	! prepare_forcing)?
466
467	! Divide shortwave into spectral bands (as in prepare_forcing)
468	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
469	CALL nemo2cice(ztmp,swvdr,'T', 1. )
470	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
471	CALL nemo2cice(ztmp,swvdf,'T', 1. )
472	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
473	CALL nemo2cice(ztmp,swidr,'T', 1. )
474	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
475	CALL nemo2cice(ztmp,swidf,'T', 1. )
476
477	ENDIF
478
479	#if defined key_asminc
480	!Ice concentration change (from assimilation)
481	ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
482	Call nemo2cice(ztmp,daice_da,'T', 1. )
483	#endif
484
485	! Snowfall
486	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
487	IF( iom_use('snowpre') ) CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit
488	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
489	CALL nemo2cice(ztmp,fsnow,'T', 1. )
490
491	! Rainfall
492	IF( iom_use('precip') ) CALL iom_put('precip', (1.0-fr_i(:,:))*(tprecip(:,:)-sprecip(:,:)) ) !!Joakim edit
493	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
494	CALL nemo2cice(ztmp,frain,'T', 1. )
495
496	! Recalculate freezing temperature and send to CICE
497	CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1))
498	CALL nemo2cice(sstfrz,Tf,'T', 1. )
499
500	! Freezing/melting potential
501	! Calculated over NEMO leapfrog timestep (hence 2*dt)
502	nfrzmlt(:,:)=rau0rcpfse3t_m(:,:)(sstfrz(:,:)-sst_m(:,:))/(2.0dt)
503	CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
504
505	! SST and SSS
506
507	CALL nemo2cice(sst_m,sst,'T', 1. )
508	CALL nemo2cice(sss_m,sss,'T', 1. )
509
510	IF( ksbc == jp_purecpl ) THEN
511	! Sea ice surface skin temperature
512	DO jl=1,ncat
513	CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
514	ENDDO
515	ENDIF
516
517	! x comp and y comp of surface ocean current
518	! U point to F point
519	DO jj=1,jpjm1
520	DO ji=1,jpi
521	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
522	ENDDO
523	ENDDO
524	CALL nemo2cice(ztmp,uocn,'F', -1. )
525
526	! V point to F point
527	DO jj=1,jpj
528	DO ji=1,jpim1
529	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
530	ENDDO
531	ENDDO
532	CALL nemo2cice(ztmp,vocn,'F', -1. )
533
534	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
535	!
536	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
537	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
538	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
539	!
540	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
541	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
542	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
543	!
544	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
545	!
546	!
547	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
548	zpice(:,:) = ssh_m(:,:)
549	ENDIF
550
551	! x comp and y comp of sea surface slope (on F points)
552	! T point to F point
553	DO jj=1,jpjm1
554	DO ji=1,jpim1
555	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj ))/e1u(ji,jj ) &
556	+ (zpice(ji+1,jj+1)-zpice(ji,jj+1))/e1u(ji,jj+1) ) &
557	* fmask(ji,jj,1)
558	ENDDO
559	ENDDO
560	CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
561
562	! T point to F point
563	DO jj=1,jpjm1
564	DO ji=1,jpim1
565	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj))/e2v(ji ,jj) &
566	+ (zpice(ji+1,jj+1)-zpice(ji+1,jj))/e2v(ji+1,jj) ) &
567	* fmask(ji,jj,1)
568	ENDDO
569	ENDDO
570	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
571
572	CALL wrk_dealloc( jpi,jpj, ztmp, zpice )
573	CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
574	!
575	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_in')
576	!
577	END SUBROUTINE cice_sbc_in
578
579
580	SUBROUTINE cice_sbc_out (kt,ksbc)
581	!!---------------------------------------------------------------------
582	!! * ROUTINE cice_sbc_out *
583	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
584	!!---------------------------------------------------------------------
585	INTEGER, INTENT( in ) :: kt ! ocean time step
586	INTEGER, INTENT( in ) :: ksbc ! surface forcing type
587
588	INTEGER :: ji, jj, jl ! dummy loop indices
589	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
590	!!---------------------------------------------------------------------
591
592	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_out')
593	!
594	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
595
596	IF( kt == nit000 ) THEN
597	IF(lwp) WRITE(numout,*)'cice_sbc_out'
598	ENDIF
599
600	! x comp of ocean-ice stress
601	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
602	ss_iou(:,:)=0.0
603	! F point to U point
604	DO jj=2,jpjm1
605	DO ji=2,jpim1
606	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
607	ENDDO
608	ENDDO
609	CALL lbc_lnk( ss_iou , 'U', -1. )
610
611	! y comp of ocean-ice stress
612	CALL cice2nemo(strocny,ztmp1,'F', -1. )
613	ss_iov(:,:)=0.0
614	! F point to V point
615
616	DO jj=1,jpjm1
617	DO ji=2,jpim1
618	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
619	ENDDO
620	ENDDO
621	CALL lbc_lnk( ss_iov , 'V', -1. )
622
623	! x and y comps of surface stress
624	! Combine wind stress and ocean-ice stress
625	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
626	! strocnx and strocny already weighted by ice fraction in CICE so not done here
627
628	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
629	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
630
631	! Also need ice/ocean stress on T points so that taum can be updated
632	! This interpolation is already done in CICE so best to use those values
633	CALL cice2nemo(strocnxT,ztmp1,'T',-1.)
634	CALL cice2nemo(strocnyT,ztmp2,'T',-1.)
635
636	! Update taum with modulus of ice-ocean stress
637	! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
638	taum(:,:)=(1.0-fr_i(:,:))taum(:,:)+fr_i(:,:)SQRT(ztmp12. + ztmp22.)
639
640	! Freshwater fluxes
641
642	IF (ksbc == jp_flx) THEN
643	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
644	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
645	! Not ideal since aice won't be the same as in the atmosphere.
646	! Better to use evap and tprecip? (but for now don't read in evap in this case)
647	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
648	ELSE IF (ksbc == jp_core) THEN
649	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
650	ELSE IF (ksbc == jp_purecpl) THEN
651	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
652	! This is currently as required with the coupling fields from the UM atmosphere
653	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
654	ENDIF
655
656	#if defined key_cice4
657	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
658	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
659	#else
660	CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
661	CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
662	#endif
663
664	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
665	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
666	! which has to be compensated for by the ocean salinity potentially going negative
667	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
668	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
669	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
670	sfx(:,:)=ztmp2(:,:)*1000.0
671	emp(:,:)=emp(:,:)-ztmp1(:,:)
672	fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
673
674	CALL lbc_lnk( emp , 'T', 1. )
675	CALL lbc_lnk( sfx , 'T', 1. )
676
677	! Solar penetrative radiation and non solar surface heat flux
678
679	! Scale qsr and qns according to ice fraction (bulk formulae only)
680
681	IF (ksbc == jp_core) THEN
682	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
683	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
684	ENDIF
685	! Take into account snow melting except for fully coupled when already in qns_tot
686	IF (ksbc == jp_purecpl) THEN
687	qsr(:,:)= qsr_tot(:,:)
688	qns(:,:)= qns_tot(:,:)
689	ELSE
690	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
691	ENDIF
692
693	! Now add in ice / snow related terms
694	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
695	#if defined key_cice4
696	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
697	#else
698	CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
699	#endif
700	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
701	CALL lbc_lnk( qsr , 'T', 1. )
702
703	DO jj=1,jpj
704	DO ji=1,jpi
705	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
706	ENDDO
707	ENDDO
708
709	#if defined key_cice4
710	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
711	#else
712	CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
713	#endif
714	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
715
716	CALL lbc_lnk( qns , 'T', 1. )
717
718	! Prepare for the following CICE time-step
719
720	CALL cice2nemo(aice,fr_i,'T', 1. )
721	IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
722	DO jl=1,ncat
723	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
724	ENDDO
725	ENDIF
726
727	! T point to U point
728	! T point to V point
729	DO jj=1,jpjm1
730	DO ji=1,jpim1
731	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
732	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
733	ENDDO
734	ENDDO
735
736	CALL lbc_lnk ( fr_iu , 'U', 1. )
737	CALL lbc_lnk ( fr_iv , 'V', 1. )
738
739	! ! embedded sea ice
740	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
741	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
742	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
743	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
744	snwice_mass_b(:,:) = snwice_mass(:,:)
745	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
746	ENDIF
747
748	#if defined key_asminc
749	! Import fresh water and salt flux due to seaice da
750	CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
751	CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
752	#endif
753
754	! Release work space
755
756	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
757	!
758	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_out')
759	!
760	END SUBROUTINE cice_sbc_out
761
762
763	SUBROUTINE cice_sbc_hadgam( kt )
764	!!---------------------------------------------------------------------
765	!! * ROUTINE cice_sbc_hadgam *
766	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
767	!!
768	!!
769	INTEGER, INTENT( in ) :: kt ! ocean time step
770	!!---------------------------------------------------------------------
771
772	INTEGER :: jl ! dummy loop index
773	INTEGER :: ierror
774
775	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_hadgam')
776	!
777	! ! =========================== !
778	! ! Prepare Coupling fields !
779	! ! =========================== !
780
781	! x and y comp of ice velocity
782
783	CALL cice2nemo(uvel,u_ice,'F', -1. )
784	CALL cice2nemo(vvel,v_ice,'F', -1. )
785
786	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
787
788	! Snow and ice thicknesses (CO_2 and CO_3)
789
790	DO jl = 1,ncat
791	CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
792	CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
793	ENDDO
794
795	#if ! defined key_cice4
796	! Meltpond fraction and depth
797	DO jl = 1,ncat
798	CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. )
799	CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. )
800	ENDDO
801	#endif
802
803
804	! If using multilayers thermodynamics in CICE then get top layer temperature
805	! and effective conductivity
806	!! When using NEMO with CICE, this change requires use of
807	!! one of the following two CICE branches:
808	!! - at CICE5.0, hadax/r1015_GSI8_with_GSI7
809	!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
810	IF (heat_capacity) THEN
811	DO jl = 1,ncat
812	CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. )
813	CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. )
814	ENDDO
815	! Convert surface temperature to Kelvin
816	tn_ice(:,:,:)=tn_ice(:,:,:)+rt0
817	ELSE
818	tn_ice(:,:,:) = 0.0
819	kn_ice(:,:,:) = 0.0
820	ENDIF
821
822	!
823	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
824	!
825	END SUBROUTINE cice_sbc_hadgam
826
827
828	SUBROUTINE cice_sbc_final
829	!!---------------------------------------------------------------------
830	!! * ROUTINE cice_sbc_final *
831	!! ** Purpose: Finalize CICE
832	!!---------------------------------------------------------------------
833
834	IF(lwp) WRITE(numout,*)'cice_sbc_final'
835
836	CALL CICE_Finalize
837
838	END SUBROUTINE cice_sbc_final
839
840	SUBROUTINE cice_sbc_force (kt)
841	!!---------------------------------------------------------------------
842	!! * ROUTINE cice_sbc_force *
843	!! ** Purpose : Provide CICE forcing from files
844	!!
845	!!---------------------------------------------------------------------
846	!! ** Method : READ monthly flux file in NetCDF files
847	!!
848	!! snowfall
849	!! rainfall
850	!! sublimation rate
851	!! topmelt (category)
852	!! botmelt (category)
853	!!
854	!! History :
855	!!----------------------------------------------------------------------
856	!! * Modules used
857	USE iom
858
859	!! * arguments
860	INTEGER, INTENT( in ) :: kt ! ocean time step
861
862	INTEGER :: ierror ! return error code
863	INTEGER :: ifpr ! dummy loop index
864	!!
865	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
866	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
867	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
868	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
869	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
870
871	!!
872	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
873	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
874	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
875	INTEGER :: ios
876	!!---------------------------------------------------------------------
877
878	! ! ====================== !
879	IF( kt == nit000 ) THEN ! First call kt=nit000 !
880	! ! ====================== !
881	! namsbc_cice is not yet in the reference namelist
882	! set file information (default values)
883	cn_dir = './' ! directory in which the model is executed
884
885	! (NB: frequency positive => hours, negative => months)
886	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
887	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
888	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
889	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
890	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
891	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
892	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
893	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
894	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
895	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
896	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
897	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
898	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
899	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
900	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
901
902	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
903	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
904	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
905
906	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
907	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
908	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
909	IF(lwm) WRITE ( numond, namsbc_cice )
910
911	! store namelist information in an array
912	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
913	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
914	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
915	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
916	slf_i(jp_bot5) = sn_bot5
917
918	! set sf structure
919	ALLOCATE( sf(jpfld), STAT=ierror )
920	IF( ierror > 0 ) THEN
921	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
922	ENDIF
923
924	DO ifpr= 1, jpfld
925	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
926	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
927	END DO
928
929	! fill sf with slf_i and control print
930	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
931	!
932	ENDIF
933
934	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
935	! ! input fields at the current time-step
936
937	! set the fluxes from read fields
938	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
939	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
940	! May be better to do this conversion somewhere else
941	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
942	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
943	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
944	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
945	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
946	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
947	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
948	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
949	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
950	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
951	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
952
953	! control print (if less than 100 time-step asked)
954	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
955	WRITE(numout,*)
956	WRITE(numout,*) ' read forcing fluxes for CICE OK'
957	CALL FLUSH(numout)
958	ENDIF
959
960	END SUBROUTINE cice_sbc_force
961
962	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
963	!!---------------------------------------------------------------------
964	!! * ROUTINE nemo2cice *
965	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
966	#if defined key_nemocice_decomp
967	!!
968	!! NEMO and CICE PE sub domains are identical, hence
969	!! there is no need to gather or scatter data from
970	!! one PE configuration to another.
971	#else
972	!! Automatically gather/scatter between
973	!! different processors and blocks
974	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
975	!! B. Gather pn into global array (png)
976	!! C. Map png into CICE global array (pcg)
977	!! D. Scatter pcg to CICE blocks (pc) + update haloes
978	#endif
979	!!---------------------------------------------------------------------
980
981	CHARACTER(len=1), INTENT( in ) :: &
982	cd_type ! nature of pn grid-point
983	! ! = T or F gridpoints
984	REAL(wp), INTENT( in ) :: &
985	psgn ! control of the sign change
986	! ! =-1 , the sign is modified following the type of b.c. used
987	! ! = 1 , no sign change
988	REAL(wp), DIMENSION(jpi,jpj) :: pn
989	#if !defined key_nemocice_decomp
990	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
991	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
992	#endif
993	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
994	INTEGER (int_kind) :: &
995	field_type, &! id for type of field (scalar, vector, angle)
996	grid_loc ! id for location on horizontal grid
997	! (center, NEcorner, Nface, Eface)
998
999	INTEGER :: ji, jj, jn ! dummy loop indices
1000
1001	! A. Ensure all haloes are filled in NEMO field (pn)
1002
1003	CALL lbc_lnk( pn , cd_type, psgn )
1004
1005	#if defined key_nemocice_decomp
1006
1007	! Copy local domain data from NEMO to CICE field
1008	pc(:,:,1)=0.0
1009	DO jj=2,ny_block-1
1010	DO ji=2,nx_block-1
1011	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
1012	ENDDO
1013	ENDDO
1014
1015	#else
1016
1017	! B. Gather pn into global array (png)
1018
1019	IF ( jpnij > 1) THEN
1020	CALL mppsync
1021	CALL mppgather (pn,0,png)
1022	CALL mppsync
1023	ELSE
1024	png(:,:,1)=pn(:,:)
1025	ENDIF
1026
1027	! C. Map png into CICE global array (pcg)
1028
1029	! Need to make sure this is robust to changes in NEMO halo rows....
1030	! (may be OK but not 100% sure)
1031
1032	IF (nproc==0) THEN
1033	! pcg(:,:)=0.0
1034	DO jn=1,jpnij
1035	DO jj=nldjt(jn),nlejt(jn)
1036	DO ji=nldit(jn),nleit(jn)
1037	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
1038	ENDDO
1039	ENDDO
1040	ENDDO
1041	DO jj=1,ny_global
1042	DO ji=1,nx_global
1043	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
1044	ENDDO
1045	ENDDO
1046	ENDIF
1047
1048	#endif
1049
1050	SELECT CASE ( cd_type )
1051	CASE ( 'T' )
1052	grid_loc=field_loc_center
1053	CASE ( 'F' )
1054	grid_loc=field_loc_NEcorner
1055	END SELECT
1056
1057	SELECT CASE ( NINT(psgn) )
1058	CASE ( -1 )
1059	field_type=field_type_vector
1060	CASE ( 1 )
1061	field_type=field_type_scalar
1062	END SELECT
1063
1064	#if defined key_nemocice_decomp
1065	! Ensure CICE halos are up to date
1066	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1067	#else
1068	! D. Scatter pcg to CICE blocks (pc) + update halos
1069	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
1070	#endif
1071
1072	END SUBROUTINE nemo2cice
1073
1074	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
1075	!!---------------------------------------------------------------------
1076	!! * ROUTINE cice2nemo *
1077	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
1078	#if defined key_nemocice_decomp
1079	!!
1080	!! NEMO and CICE PE sub domains are identical, hence
1081	!! there is no need to gather or scatter data from
1082	!! one PE configuration to another.
1083	#else
1084	!! Automatically deal with scatter/gather between
1085	!! different processors and blocks
1086	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
1087	!! B. Map pcg into NEMO global array (png)
1088	!! C. Scatter png into NEMO field (pn) for each processor
1089	!! D. Ensure all haloes are filled in pn
1090	#endif
1091	!!---------------------------------------------------------------------
1092
1093	CHARACTER(len=1), INTENT( in ) :: &
1094	cd_type ! nature of pn grid-point
1095	! ! = T or F gridpoints
1096	REAL(wp), INTENT( in ) :: &
1097	psgn ! control of the sign change
1098	! ! =-1 , the sign is modified following the type of b.c. used
1099	! ! = 1 , no sign change
1100	REAL(wp), DIMENSION(jpi,jpj) :: pn
1101
1102	#if defined key_nemocice_decomp
1103	INTEGER (int_kind) :: &
1104	field_type, & ! id for type of field (scalar, vector, angle)
1105	grid_loc ! id for location on horizontal grid
1106	! (center, NEcorner, Nface, Eface)
1107	#else
1108	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
1109	#endif
1110
1111	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
1112
1113	INTEGER :: ji, jj, jn ! dummy loop indices
1114
1115
1116	#if defined key_nemocice_decomp
1117
1118	SELECT CASE ( cd_type )
1119	CASE ( 'T' )
1120	grid_loc=field_loc_center
1121	CASE ( 'F' )
1122	grid_loc=field_loc_NEcorner
1123	END SELECT
1124
1125	SELECT CASE ( NINT(psgn) )
1126	CASE ( -1 )
1127	field_type=field_type_vector
1128	CASE ( 1 )
1129	field_type=field_type_scalar
1130	END SELECT
1131
1132	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1133
1134
1135	pn(:,:)=0.0
1136	DO jj=1,jpjm1
1137	DO ji=1,jpim1
1138	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
1139	ENDDO
1140	ENDDO
1141
1142	#else
1143
1144	! A. Gather CICE blocks (pc) into global array (pcg)
1145
1146	CALL gather_global(pcg, pc, 0, distrb_info)
1147
1148	! B. Map pcg into NEMO global array (png)
1149
1150	! Need to make sure this is robust to changes in NEMO halo rows....
1151	! (may be OK but not spent much time thinking about it)
1152	! Note that non-existent pcg elements may be used below, but
1153	! the lbclnk call on pn will replace these with sensible values
1154
1155	IF (nproc==0) THEN
1156	png(:,:,:)=0.0
1157	DO jn=1,jpnij
1158	DO jj=nldjt(jn),nlejt(jn)
1159	DO ji=nldit(jn),nleit(jn)
1160	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
1161	ENDDO
1162	ENDDO
1163	ENDDO
1164	ENDIF
1165
1166	! C. Scatter png into NEMO field (pn) for each processor
1167
1168	IF ( jpnij > 1) THEN
1169	CALL mppsync
1170	CALL mppscatter (png,0,pn)
1171	CALL mppsync
1172	ELSE
1173	pn(:,:)=png(:,:,1)
1174	ENDIF
1175
1176	#endif
1177
1178	! D. Ensure all haloes are filled in pn
1179
1180	CALL lbc_lnk( pn , cd_type, psgn )
1181
1182	END SUBROUTINE cice2nemo
1183
1184	#else
1185	!!----------------------------------------------------------------------
1186	!! Default option Dummy module NO CICE sea-ice model
1187	!!----------------------------------------------------------------------
1188	!! $Id$
1189	CONTAINS
1190
1191	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
1192	IMPLICIT NONE
1193	INTEGER, INTENT(in) :: kt ! ocean time-step index
1194	INTEGER, INTENT(in) :: ksbc ! surface forcing type
1195	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
1196	END SUBROUTINE sbc_ice_cice
1197
1198	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
1199	IMPLICIT NONE
1200	INTEGER, INTENT(in) :: ksbc ! surface forcing type
1201	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
1202	END SUBROUTINE cice_sbc_init
1203
1204	SUBROUTINE cice_sbc_final ! Dummy routine
1205	IMPLICIT NONE
1206	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
1207	END SUBROUTINE cice_sbc_final
1208
1209	#endif
1210
1211	!!======================================================================
1212	END MODULE sbcice_cice

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

Download in other formats: