Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcice_cice.F90 @ 14627

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

Download in other formats: