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

Context Navigation

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5797

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

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

Download in other formats: