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_r5518_GO6_package_05sprecip/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_05sprecip/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9281

Last change on this file since 9281 was 9281, checked in by anaguiar, 6 years ago
Correct default value and add print out the value of snow precip
File size: 47.6 KB

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

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

Download in other formats: