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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8819

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

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

Download in other formats: