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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_w_keyasminc/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 7749

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

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

Download in other formats: