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

Context Navigation

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

Last change on this file since 8384 was 8384, checked in by jamrae, 7 years ago
Implemented changes to pass sea ice form drag from CICE to coupler.
File size: 46.9 KB

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

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

Download in other formats: