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 @ 5393

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

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

Download in other formats: