Context Navigation

sbcice_cice.F90 @ 11107

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

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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11107

Download in other formats: