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

Context Navigation

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9695

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

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

Download in other formats: