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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8190

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

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

Download in other formats: