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

Context Navigation

source: branches/UKMO/dev_r5518_CICE_coupling_GSI7/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5612

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

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

Download in other formats: