Context Navigation

sbcice_cice.F90 @ 9383

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

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

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

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9383

Download in other formats: