Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

sbcice_cice.F90

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

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

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

Context Navigation

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

Download in other formats: