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

Context Navigation

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5796

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

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

Download in other formats: