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

Context Navigation

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_with_key_asminc/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6172

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

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

Download in other formats: