Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcice_cice.F90 @ 11738

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

Download in other formats: