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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcice_cice.F90 in branches/UKMO/dev_r5518_GO6_package_r8638_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

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

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

Line
1	MODULE sbcice_cice
2	!!======================================================================
3	!! * MODULE sbcice_cice *
4	!! To couple with sea ice model CICE (LANL)
5	!!=====================================================================
6	#if defined key_cice
7	!!----------------------------------------------------------------------
8	!! 'key_cice' : CICE sea-ice model
9	!!----------------------------------------------------------------------
10	!! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area
11	!!
12	!!
13	!!----------------------------------------------------------------------
14	USE oce ! ocean dynamics and tracers
15	USE dom_oce ! ocean space and time domain
16	USE domvvl
17	USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
18	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
19	USE in_out_manager ! I/O manager
20	USE iom, ONLY : iom_put,iom_use ! I/O manager library !!Joakim edit
21	USE lib_mpp ! distributed memory computing library
22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
23	USE wrk_nemo ! work arrays
24	USE timing ! Timing
25	USE daymod ! calendar
26	USE fldread ! read input fields
27	USE sbc_oce ! Surface boundary condition: ocean fields
28	USE sbc_ice ! Surface boundary condition: ice fields
29	USE sbcblk_core ! Surface boundary condition: CORE bulk
30	USE sbccpl
31
32	USE ice_kinds_mod
33	USE ice_blocks
34	USE ice_domain
35	USE ice_domain_size
36	USE ice_boundary
37	USE ice_constants
38	USE ice_gather_scatter
39	USE ice_calendar, only: dt
40	# if defined key_cice4
41	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
42	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
43	strocnxT,strocnyT, &
44	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
45	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
46	flatn_f,fsurfn_f,fcondtopn_f, &
47	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
48	swvdr,swvdf,swidr,swidf,Tf
49	USE ice_therm_vertical, only: calc_Tsfc
50	#else
51	USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
52	vsnon,vice,vicen,nt_Tsfc
53	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
54	strocnxT,strocnyT, &
55	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai, &
56	fresh_ai,fhocn_ai,fswthru_ai,frzmlt, &
57	flatn_f,fsurfn_f,fcondtopn_f, &
58	#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	CALL cice2nemo(z0atmskin,rough_ice_skin,'T', 1. )
819
820	!
821	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
822	!
823	END SUBROUTINE cice_sbc_hadgam
824
825
826	SUBROUTINE cice_sbc_final
827	!!---------------------------------------------------------------------
828	!! * ROUTINE cice_sbc_final *
829	!! ** Purpose: Finalize CICE
830	!!---------------------------------------------------------------------
831
832	IF(lwp) WRITE(numout,*)'cice_sbc_final'
833
834	CALL CICE_Finalize
835
836	END SUBROUTINE cice_sbc_final
837
838	SUBROUTINE cice_sbc_force (kt)
839	!!---------------------------------------------------------------------
840	!! * ROUTINE cice_sbc_force *
841	!! ** Purpose : Provide CICE forcing from files
842	!!
843	!!---------------------------------------------------------------------
844	!! ** Method : READ monthly flux file in NetCDF files
845	!!
846	!! snowfall
847	!! rainfall
848	!! sublimation rate
849	!! topmelt (category)
850	!! botmelt (category)
851	!!
852	!! History :
853	!!----------------------------------------------------------------------
854	!! * Modules used
855	USE iom
856
857	!! * arguments
858	INTEGER, INTENT( in ) :: kt ! ocean time step
859
860	INTEGER :: ierror ! return error code
861	INTEGER :: ifpr ! dummy loop index
862	!!
863	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
864	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
865	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
866	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
867	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
868
869	!!
870	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
871	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
872	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
873	INTEGER :: ios
874	!!---------------------------------------------------------------------
875
876	! ! ====================== !
877	IF( kt == nit000 ) THEN ! First call kt=nit000 !
878	! ! ====================== !
879	! namsbc_cice is not yet in the reference namelist
880	! set file information (default values)
881	cn_dir = './' ! directory in which the model is executed
882
883	! (NB: frequency positive => hours, negative => months)
884	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! landmask
885	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! file
886	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' , '' )
887	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' , '' )
888	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' , '' )
889	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
890	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
891	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
892	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
893	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
894	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' , '' )
895	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' , '' )
896	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' , '' )
897	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' , '' )
898	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' , '' )
899
900	REWIND( numnam_ref ) ! Namelist namsbc_cice in reference namelist :
901	READ ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
902	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in reference namelist', lwp )
903
904	REWIND( numnam_cfg ) ! Namelist namsbc_cice in configuration namelist : Parameters of the run
905	READ ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
906	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist', lwp )
907	IF(lwm) WRITE ( numond, namsbc_cice )
908
909	! store namelist information in an array
910	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
911	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
912	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
913	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
914	slf_i(jp_bot5) = sn_bot5
915
916	! set sf structure
917	ALLOCATE( sf(jpfld), STAT=ierror )
918	IF( ierror > 0 ) THEN
919	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
920	ENDIF
921
922	DO ifpr= 1, jpfld
923	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
924	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
925	END DO
926
927	! fill sf with slf_i and control print
928	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
929	!
930	ENDIF
931
932	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
933	! ! input fields at the current time-step
934
935	! set the fluxes from read fields
936	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
937	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
938	! May be better to do this conversion somewhere else
939	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
940	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
941	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
942	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
943	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
944	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
945	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
946	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
947	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
948	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
949	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
950
951	! control print (if less than 100 time-step asked)
952	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
953	WRITE(numout,*)
954	WRITE(numout,*) ' read forcing fluxes for CICE OK'
955	CALL FLUSH(numout)
956	ENDIF
957
958	END SUBROUTINE cice_sbc_force
959
960	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
961	!!---------------------------------------------------------------------
962	!! * ROUTINE nemo2cice *
963	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
964	#if defined key_nemocice_decomp
965	!!
966	!! NEMO and CICE PE sub domains are identical, hence
967	!! there is no need to gather or scatter data from
968	!! one PE configuration to another.
969	#else
970	!! Automatically gather/scatter between
971	!! different processors and blocks
972	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
973	!! B. Gather pn into global array (png)
974	!! C. Map png into CICE global array (pcg)
975	!! D. Scatter pcg to CICE blocks (pc) + update haloes
976	#endif
977	!!---------------------------------------------------------------------
978
979	CHARACTER(len=1), INTENT( in ) :: &
980	cd_type ! nature of pn grid-point
981	! ! = T or F gridpoints
982	REAL(wp), INTENT( in ) :: &
983	psgn ! control of the sign change
984	! ! =-1 , the sign is modified following the type of b.c. used
985	! ! = 1 , no sign change
986	REAL(wp), DIMENSION(jpi,jpj) :: pn
987	#if !defined key_nemocice_decomp
988	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
989	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
990	#endif
991	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
992	INTEGER (int_kind) :: &
993	field_type, &! id for type of field (scalar, vector, angle)
994	grid_loc ! id for location on horizontal grid
995	! (center, NEcorner, Nface, Eface)
996
997	INTEGER :: ji, jj, jn ! dummy loop indices
998
999	! A. Ensure all haloes are filled in NEMO field (pn)
1000
1001	CALL lbc_lnk( pn , cd_type, psgn )
1002
1003	#if defined key_nemocice_decomp
1004
1005	! Copy local domain data from NEMO to CICE field
1006	pc(:,:,1)=0.0
1007	DO jj=2,ny_block-1
1008	DO ji=2,nx_block-1
1009	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
1010	ENDDO
1011	ENDDO
1012
1013	#else
1014
1015	! B. Gather pn into global array (png)
1016
1017	IF ( jpnij > 1) THEN
1018	CALL mppsync
1019	CALL mppgather (pn,0,png)
1020	CALL mppsync
1021	ELSE
1022	png(:,:,1)=pn(:,:)
1023	ENDIF
1024
1025	! C. Map png into CICE global array (pcg)
1026
1027	! Need to make sure this is robust to changes in NEMO halo rows....
1028	! (may be OK but not 100% sure)
1029
1030	IF (nproc==0) THEN
1031	! pcg(:,:)=0.0
1032	DO jn=1,jpnij
1033	DO jj=nldjt(jn),nlejt(jn)
1034	DO ji=nldit(jn),nleit(jn)
1035	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
1036	ENDDO
1037	ENDDO
1038	ENDDO
1039	DO jj=1,ny_global
1040	DO ji=1,nx_global
1041	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
1042	ENDDO
1043	ENDDO
1044	ENDIF
1045
1046	#endif
1047
1048	SELECT CASE ( cd_type )
1049	CASE ( 'T' )
1050	grid_loc=field_loc_center
1051	CASE ( 'F' )
1052	grid_loc=field_loc_NEcorner
1053	END SELECT
1054
1055	SELECT CASE ( NINT(psgn) )
1056	CASE ( -1 )
1057	field_type=field_type_vector
1058	CASE ( 1 )
1059	field_type=field_type_scalar
1060	END SELECT
1061
1062	#if defined key_nemocice_decomp
1063	! Ensure CICE halos are up to date
1064	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1065	#else
1066	! D. Scatter pcg to CICE blocks (pc) + update halos
1067	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
1068	#endif
1069
1070	END SUBROUTINE nemo2cice
1071
1072	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
1073	!!---------------------------------------------------------------------
1074	!! * ROUTINE cice2nemo *
1075	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
1076	#if defined key_nemocice_decomp
1077	!!
1078	!! NEMO and CICE PE sub domains are identical, hence
1079	!! there is no need to gather or scatter data from
1080	!! one PE configuration to another.
1081	#else
1082	!! Automatically deal with scatter/gather between
1083	!! different processors and blocks
1084	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
1085	!! B. Map pcg into NEMO global array (png)
1086	!! C. Scatter png into NEMO field (pn) for each processor
1087	!! D. Ensure all haloes are filled in pn
1088	#endif
1089	!!---------------------------------------------------------------------
1090
1091	CHARACTER(len=1), INTENT( in ) :: &
1092	cd_type ! nature of pn grid-point
1093	! ! = T or F gridpoints
1094	REAL(wp), INTENT( in ) :: &
1095	psgn ! control of the sign change
1096	! ! =-1 , the sign is modified following the type of b.c. used
1097	! ! = 1 , no sign change
1098	REAL(wp), DIMENSION(jpi,jpj) :: pn
1099
1100	#if defined key_nemocice_decomp
1101	INTEGER (int_kind) :: &
1102	field_type, & ! id for type of field (scalar, vector, angle)
1103	grid_loc ! id for location on horizontal grid
1104	! (center, NEcorner, Nface, Eface)
1105	#else
1106	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
1107	#endif
1108
1109	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
1110
1111	INTEGER :: ji, jj, jn ! dummy loop indices
1112
1113
1114	#if defined key_nemocice_decomp
1115
1116	SELECT CASE ( cd_type )
1117	CASE ( 'T' )
1118	grid_loc=field_loc_center
1119	CASE ( 'F' )
1120	grid_loc=field_loc_NEcorner
1121	END SELECT
1122
1123	SELECT CASE ( NINT(psgn) )
1124	CASE ( -1 )
1125	field_type=field_type_vector
1126	CASE ( 1 )
1127	field_type=field_type_scalar
1128	END SELECT
1129
1130	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1131
1132
1133	pn(:,:)=0.0
1134	DO jj=1,jpjm1
1135	DO ji=1,jpim1
1136	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
1137	ENDDO
1138	ENDDO
1139
1140	#else
1141
1142	! A. Gather CICE blocks (pc) into global array (pcg)
1143
1144	CALL gather_global(pcg, pc, 0, distrb_info)
1145
1146	! B. Map pcg into NEMO global array (png)
1147
1148	! Need to make sure this is robust to changes in NEMO halo rows....
1149	! (may be OK but not spent much time thinking about it)
1150	! Note that non-existent pcg elements may be used below, but
1151	! the lbclnk call on pn will replace these with sensible values
1152
1153	IF (nproc==0) THEN
1154	png(:,:,:)=0.0
1155	DO jn=1,jpnij
1156	DO jj=nldjt(jn),nlejt(jn)
1157	DO ji=nldit(jn),nleit(jn)
1158	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
1159	ENDDO
1160	ENDDO
1161	ENDDO
1162	ENDIF
1163
1164	! C. Scatter png into NEMO field (pn) for each processor
1165
1166	IF ( jpnij > 1) THEN
1167	CALL mppsync
1168	CALL mppscatter (png,0,pn)
1169	CALL mppsync
1170	ELSE
1171	pn(:,:)=png(:,:,1)
1172	ENDIF
1173
1174	#endif
1175
1176	! D. Ensure all haloes are filled in pn
1177
1178	CALL lbc_lnk( pn , cd_type, psgn )
1179
1180	END SUBROUTINE cice2nemo
1181
1182	#else
1183	!!----------------------------------------------------------------------
1184	!! Default option Dummy module NO CICE sea-ice model
1185	!!----------------------------------------------------------------------
1186	!! $Id$
1187	CONTAINS
1188
1189	SUBROUTINE sbc_ice_cice ( kt, ksbc ) ! Dummy routine
1190	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
1191	END SUBROUTINE sbc_ice_cice
1192
1193	SUBROUTINE cice_sbc_init (ksbc) ! Dummy routine
1194	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
1195	END SUBROUTINE cice_sbc_init
1196
1197	SUBROUTINE cice_sbc_final ! Dummy routine
1198	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
1199	END SUBROUTINE cice_sbc_final
1200
1201	#endif
1202
1203	!!======================================================================
1204	END MODULE sbcice_cice

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

Download in other formats: