New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcice_cice.F90 in branches/UKMO/dev_r5377_UKMO1_CICE_coupling_GSI7/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5377_UKMO1_CICE_coupling_GSI7/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5522

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

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

Download in other formats: