Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

sbcice_cice.F90

Last change on this file was 12613, checked in by dcarneir, 4 years ago
Including sea ice thickness in GO6 branch
File size: 49.2 KB

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

Download in other formats: