Context Navigation

sbcice_cice.F90 @ 9717

Last change on this file since 9717 was 9717, checked in by jamrae, 6 years ago
Made changes for sea ice form drag coupling.
File size: 47.7 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_r8638_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9717

Download in other formats: