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

Context Navigation

source: branches/2015/dev_r5177_CNRS4_stopar/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5304

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

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

Download in other formats: