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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10770

Last change on this file since 10770 was 10770, checked in by andmirek, 6 years ago
GMED 450 changes for writing to numout
File size: 47.8 KB

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

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

Download in other formats: