Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

sbcice_cice.F90

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

Download in other formats: