Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcice_cice.F90 @ 8975

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

Download in other formats: